您的位置:首页 > 运维架构

灰度共生矩阵及相关特征值的计算——opencv

2017-06-19 11:13 525 查看
#include<iostream>
#include<opencv2/highgui.hpp>
#include<opencv2/core.hpp>
#include<opencv2/imgcodecs.hpp>
#include<opencv2/opencv.hpp>
using namespace std;
using namespace cv;

const int gray_level = 16;

void getglcm_horison(Mat& input, Mat& dst)//0度灰度共生矩阵
{
Mat src=input;
CV_Assert(1 == src.channels());
src.convertTo(src, CV_32S);
int height = src.rows;
int width = src.cols;
int max_gray_level=0;
for (int j = 0; j < height; j++)//寻找像素灰度最大值
{
int* srcdata = src.ptr<int>(j);
for (int i = 0; i < width; i++)
{
if (srcdata[i] > max_gray_level)
{
max_gray_level = srcdata[i];
}
}
}
max_gray_level++;//像素灰度最大值加1即为该矩阵所拥有的灰度级数
if (max_gray_level > 16)//若灰度级数大于16,则将图像的灰度级缩小至16级,减小灰度共生矩阵的大小。
{
for (int i = 0; i < height; i++)
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width; j++)
{
srcdata[j] = (int)srcdata[j] / gray_level;
}
}

dst.create(gray_level, gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height; i++)
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width - 1; j++)
{
int rows = srcdata[j];
int cols = srcdata[j + 1];
dst.ptr<int>(rows)[cols]++;
}
}
}
else//若灰度级数小于16,则生成相应的灰度共生矩阵
{
dst.create(max_gray_level, max_gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height; i++)
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width - 1; j++)
{
int rows = srcdata[j];
int cols = srcdata[j + 1];
dst.ptr<int>(rows)[cols]++;
}
}
}
}

void getglcm_vertical(Mat& input, Mat& dst)//90度灰度共生矩阵
{
Mat src = input;
CV_Assert(1 == src.channels());
src.convertTo(src, CV_32S);
int height = src.rows;
int width = src.cols;
int max_gray_level = 0;
for (int j = 0; j < height; j++)
{
int* srcdata = src.ptr<int>(j);
for (int i = 0; i < width; i++)
{
if (srcdata[i] > max_gray_level)
{
max_gray_level = srcdata[i];
}
}
}
max_gray_level++;
if (max_gray_level > 16)
{
for (int i = 0; i < height; i++)//将图像的灰度级缩小至16级,减小灰度共生矩阵的大小。
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width; j++)
{
srcdata[j] = (int)srcdata[j] / gray_level;
}
}

dst.create(gray_level, gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height-1; i++)
{
int*srcdata = src.ptr<int>(i);
int*srcdata1 = src.ptr<int>(i+1);
for (int j = 0; j < width ; j++)
{
int rows = srcdata[j];
int cols = srcdata1[j];
dst.ptr<int>(rows)[cols]++;
}
}
}
else
{
dst.create(max_gray_level, max_gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height-1; i++)
{
int*srcdata = src.ptr<int>(i);
int*srcdata1 = src.ptr<int>(i + 1);
for (int j = 0; j < width; j++)
{
int rows = srcdata[j];
int cols = srcdata1[j];
dst.ptr<int>(rows)[cols]++;
}
}
}
}

void getglcm_45(Mat& input, Mat& dst)//45度灰度共生矩阵
{
Mat src = input;
CV_Assert(1 == src.channels());
src.convertTo(src, CV_32S);
int height = src.rows;
int width = src.cols;
int max_gray_level = 0;
for (int j = 0; j < height; j++)
{
int* srcdata = src.ptr<int>(j);
for (int i = 0; i < width; i++)
{
if (srcdata[i] > max_gray_level)
{
max_gray_level = srcdata[i];
}
}
}
max_gray_level++;
if (max_gray_level > 16)
{
for (int i = 0; i < height; i++)//将图像的灰度级缩小至16级,减小灰度共生矩阵的大小。
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width; j++)
{
srcdata[j] = (int)srcdata[j] / gray_level;
}
}

dst.create(gray_level, gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height - 1; i++)
{
int*srcdata = src.ptr<int>(i);
int*srcdata1 = src.ptr<int>(i + 1);
for (int j = 0; j < width-1; j++)
{
int rows = srcdata[j];
int cols = srcdata1[j+1];
dst.ptr<int>(rows)[cols]++;
}
}
}
else
{
dst.create(max_gray_level, max_gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height - 1; i++)
{
int*srcdata = src.ptr<int>(i);
int*srcdata1 = src.ptr<int>(i + 1);
for (int j = 0; j < width-1; j++)
{
int rows = srcdata[j];
int cols = srcdata1[j+1];
dst.ptr<int>(rows)[cols]++;
}
}
}
}

void getglcm_135(Mat& input, Mat& dst)//135度灰度共生矩阵
{
Mat src = input;
CV_Assert(1 == src.channels());
src.convertTo(src, CV_32S);
int height = src.rows;
int width = src.cols;
int max_gray_level = 0;
for (int j = 0; j < height; j++)
{
int* srcdata = src.ptr<int>(j);
for (int i = 0; i < width; i++)
{
if (srcdata[i] > max_gray_level)
{
max_gray_level = srcdata[i];
}
}
}
max_gray_level++;
if (max_gray_level > 16)
{
for (int i = 0; i < height; i++)//将图像的灰度级缩小至16级,减小灰度共生矩阵的大小。
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width; j++)
{
srcdata[j] = (int)srcdata[j] / gray_level;
}
}

dst.create(gray_level, gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height - 1; i++)
{
int*srcdata = src.ptr<int>(i);
int*srcdata1 = src.ptr<int>(i + 1);
for (int j = 1; j < width; j++)
{
int rows = srcdata[j];
int cols = srcdata1[j-1];
dst.ptr<int>(rows)[cols]++;
}
}
}
else
{
dst.create(max_gray_level, max_gray_level, CV_32SC1);
dst = Scalar::all(0);
for (int i = 0; i < height - 1; i++)
{
int*srcdata = src.ptr<int>(i);
int*srcdata1 = src.ptr<int>(i + 1);
for (int j = 1; j < width; j++)
{
int rows = srcdata[j];
int cols = srcdata1[j-1];
dst.ptr<int>(rows)[cols]++;
}
}
}
}

void feature_computer(Mat&src, double& Asm, double& Eng, double& Con, double& Idm)//计算特征值
{
int height = src.rows;
int width = src.cols;
int total = 0;
for (int i = 0; i < height; i++)
{
int*srcdata = src.ptr<int>(i);
for (int j = 0; j < width; j++)
{
total += srcdata[j];//求图像所有像素的灰度值的和
}
}

Mat copy;
copy.create(height, width, CV_64FC1);
for (int i = 0; i < height; i++)
{
int*srcdata = src.ptr<int>(i);
double*copydata = copy.ptr<double>(i);
for (int j = 0; j < width; j++)
{
copydata[j]=(double)srcdata[j]/(double)total;//图像每一个像素的的值除以像素总和
}
}

for (int i = 0; i < height; i++)
{
double*srcdata = copy.ptr<double>(i);
for (int j = 0; j < width; j++)
{
Asm += srcdata[j] * srcdata[j];//能量
if (srcdata[j]>0)
Eng -= srcdata[j] * log(srcdata[j]);//熵
Con += (double)(i - j)*(double)(i - j)*srcdata[j];//对比度
Idm += srcdata[j] / (1 + (double)(i - j)*(double)(i - j));//逆差矩
}
}
}

int main()
{
Mat dst_horison, dst_vertical, dst_45, dst_135;

Mat src = imread("C:\\Users\\aoe\\Desktop\\Visual C+\\Visual C+\\chapter08\\pic\\healthy\\201.bmp");
if (src.empty())
{
return -1;
}
Mat src_gray;
//src.create(src.size(), CV_8UC1);
//src_gray = Scalar::all(0);
cvtColor(src, src_gray, COLOR_BGR2GRAY);

//src =( Mat_<int>(6, 6) << 0, 1, 2, 3, 0, 1, 1, 2, 3, 0, 1, 2, 2, 3, 0, 1, 2, 3, 3, 0, 1, 2, 3, 0, 0, 1, 2, 3, 0, 1, 1, 2, 3, 0, 1, 2 );
//src = (Mat_<int>(4, 4) << 1, 17, 0, 3,3,2,20,5,26,50,1,2,81,9,25,1);
getglcm_horison(src_gray, dst_horison);
getglcm_vertical(src_gray, dst_vertical);
getglcm_45(src_gray, dst_45);
getglcm_135(src_gray, dst_135);

double eng_horison=0, con_horison=0, idm_horison=0, asm_horison=0;
feature_computer(dst_horison, asm_horison, eng_horison, con_horison, idm_horison);

cout << "asm_horison:" << asm_horison << endl;
cout << "eng_horison:" << eng_horison << endl;
cout << "con_horison:" << con_horison << endl;
cout << "idm_horison:" << idm_horison << endl;

/*  for (int i = 0; i < dst_horison.rows; i++)
{
int *data = dst_horison.ptr<int>(i);
for (int j = 0; j < dst_horison.cols; j++)
{
cout << data[j] << " ";
}
cout << endl;
}
cout << endl;

for (int i = 0; i < dst_vertical.rows; i++)
{
int *data = dst_vertical.ptr<int>(i);
for (int j = 0; j < dst_vertical.cols; j++)
{
cout << data[j] << " ";
}
cout << endl;
}
cout << endl;

for (int i = 0; i < dst_45.rows; i++)
{
int *data = dst_45.ptr<int>(i);
for (int j = 0; j < dst_45.cols; j++)
{
cout << data[j] << " ";
}
cout << endl;
}

cout << endl;

for (int i = 0; i < dst_135.rows; i++)
{
int *data = dst_135.ptr<int>(i);
for (int j = 0; j < dst_135.cols; j++)
{
cout << data[j] << " ";
}
cout << endl;
}*/
system("pause");
return 0;
}

参考:

灰度共生矩阵的定义与理解:http://www.cnblogs.com/xiangshancuizhu/archive/2011/07/24/2115266.html

OpenCV实现:

http://blog.csdn.net/cxf7394373/article/details/6988229

http://download.csdn.net/download/sxnzxz/3419181

灰度共生矩阵
        灰度共生矩阵定义为像素对的联合分布概率,是一个对称矩阵,它不仅反映图像灰度在相邻的方向、相邻间隔、变化幅度的综合信息,但也反映了相同的灰度级像素之间的位置分布特征,是计算纹理特征的基础。

       设f(x,y)为一幅数字图像,其大小为M×N,灰度级别为Ng,则满足一定空间关系的灰度共生矩阵为:



        其中#(x)表示集合x中的元素个数,显然P为Ng×Ng的矩阵,若(x1,y1)与(x2,y2)间距离为d,两者与坐标横轴的夹角为θ,则可以得到各种间距及角度的灰度共生矩阵(i,j,d,θ)。其中元素(i,j)的值表示一个灰度为i,另一个灰度为j的两个相距为d的像素对在角的方向上出现的次数。

在计算得到共生矩阵之后,往往不是直接应用计算的灰度共生矩阵,而是在此基础上计算纹理特征量,我们经常用反差、能量、熵、相关性等特征量来表示纹理特征。
        (1) 反差:又称为对比度,度量矩阵的值是如何分布和图像中局部变化的多少,反应了图像的清晰度和纹理的沟纹深浅。纹理的沟纹越深,反差越大,效果清晰;反之,对比值小,则沟纹浅,效果模糊。 



        (2) 能量:是灰度共生矩阵各元素值的平方和,是对图像纹理的灰度变化稳定程度的度量,反应了图像灰度分布均匀程度和纹理粗细度。能量值大表明当前纹理是一种规则变化较为稳定的纹理。          



        (3) 熵:是图像包含信息量的随机性度量。当共生矩阵中所有值均相等或者像素值表现出最大的随机性时,熵最大;因此熵值表明了图像灰度分布的复杂程度,熵值越大,图像越复杂。             



        (4) 相关性:也称为同质性,用来度量图像的灰度级在行或列方向上的相似程度,因此值的大小反应了局部灰度相关性,值越大,相关性也越大。



应用

    由上面的叙述知道,可以根据各种间距和角度计算灰度共生矩阵,下面程序中给定了间距,根据传入的参数计算:

[cpp]
view plain
copy

#define GLCM_DIS 3  //灰度共生矩阵的统计距离  
#define GLCM_CLASS 16 //计算灰度共生矩阵的图像灰度值等级化  
#define GLCM_ANGLE_HORIZATION 0  //水平  
#define GLCM_ANGLE_VERTICAL   1  //垂直  
#define GLCM_ANGLE_DIGONAL    2  //对角  
int calGLCM(IplImage* bWavelet,int angleDirection,double* featureVector)  
{  
    int i,j;  
    int width,height;  
  
    if(NULL == bWavelet)  
        return 1;  
  
    width = bWavelet->width;  
    height = bWavelet->height;  
  
    int * glcm = new int[GLCM_CLASS * GLCM_CLASS];  
    int * histImage = new int[width * height];  
  
    if(NULL == glcm || NULL == histImage)  
        return 2;  
  
    //灰度等级化---分GLCM_CLASS个等级  
    uchar *data =(uchar*) bWavelet->imageData;  
    for(i = 0;i < height;i++){  
        for(j = 0;j < width;j++){  
            histImage[i * width + j] = (int)(data[bWavelet->widthStep * i + j] * GLCM_CLASS / 256);  
        }  
    }  
  
    //初始化共生矩阵  
    for (i = 0;i < GLCM_CLASS;i++)  
        for (j = 0;j < GLCM_CLASS;j++)  
            glcm[i * GLCM_CLASS + j] = 0;  
  
    //计算灰度共生矩阵  
    int w,k,l;  
    //水平方向  
    if(angleDirection == GLCM_ANGLE_HORIZATION)  
    {  
        for (i = 0;i < height;i++)  
        {  
            for (j = 0;j < width;j++)  
            {  
                l = histImage[i * width + j];  
                if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width)  
                {  
                    k = histImage[i * width + j + GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
                if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width)  
                {  
                    k = histImage[i * width + j - GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
            }  
        }  
    }  
    //垂直方向  
    else if(angleDirection == GLCM_ANGLE_VERTICAL)  
    {  
        for (i = 0;i < height;i++)  
        {  
            for (j = 0;j < width;j++)  
            {  
                l = histImage[i * width + j];  
                if(i + GLCM_DIS >= 0 && i + GLCM_DIS < height)   
                {  
                    k = histImage[(i + GLCM_DIS) * width + j];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
                if(i - GLCM_DIS >= 0 && i - GLCM_DIS < height)   
                {  
                    k = histImage[(i - GLCM_DIS) * width + j];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
            }  
        }  
    }  
    //对角方向  
    else if(angleDirection == GLCM_ANGLE_DIGONAL)  
    {  
        for (i = 0;i < height;i++)  
        {  
            for (j = 0;j < width;j++)  
            {  
                l = histImage[i * width + j];  
  
                if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width && i + GLCM_DIS >= 0 && i + GLCM_DIS < height)  
                {  
                    k = histImage[(i + GLCM_DIS) * width + j + GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
                if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width && i - GLCM_DIS >= 0 && i - GLCM_DIS < height)  
                {  
                    k = histImage[(i - GLCM_DIS) * width + j - GLCM_DIS];  
                    glcm[l * GLCM_CLASS + k]++;  
                }  
            }  
        }  
    }  
  
    //计算特征值  
    double entropy = 0,energy = 0,contrast = 0,homogenity = 0;  
    for (i = 0;i < GLCM_CLASS;i++)  
    {  
        for (j = 0;j < GLCM_CLASS;j++)  
        {  
            //熵  
            if(glcm[i * GLCM_CLASS + j] > 0)  
                entropy -= glcm[i * GLCM_CLASS + j] * log10(double(glcm[i * GLCM_CLASS + j]));  
            //能量  
            energy += glcm[i * GLCM_CLASS + j] * glcm[i * GLCM_CLASS + j];  
            //对比度  
            contrast += (i - j) * (i - j) * glcm[i * GLCM_CLASS + j];  
            //一致性  
            homogenity += 1.0 / (1 + (i - j) * (i - j)) * glcm[i * GLCM_CLASS + j];  
        }  
    }  
    //返回特征值  
    i = 0;  
    featureVector[i++] = entropy;  
    featureVector[i++] = energy;  
    featureVector[i++] = contrast;  
    featureVector[i++] = homogenity;  
  
    delete[] glcm;  
    delete[] histImage;  
    return 0;  


cvGetGLCMDescriptorStatistics
http://bbs.csdn.net/topics/360141907
原文地址:http://blog.csdn.net/yanxiaopan/article/details/52356777
http://blog.csdn.net/cxf7394373/article/details/6988229
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: