利用灰度共生矩陣提取影象紋理特徵
阿新 • • 發佈：2019-01-07
1. 灰度共生矩陣概念
灰度共生矩陣定義為畫素對的聯合概率分佈，是一個對稱矩陣，它不僅反映影象灰度在相鄰的方向、相鄰間隔、變化幅度的綜合資訊，也反映了相同的灰度級畫素之間的位置分佈特徵，是計算紋理特徵的基礎。
在影象中任意取一點(x,y)及偏離它的一點(x+a,y+b)(其中，a、b為整數，人為定義)構成點對。設該點對的灰度值為(f1,f2)，再令點(x,y)在整幅影象上移動，則會得到不同的(f1,f2)值。
設影象的最大灰度級為L，則f1與f2的組合共有L*L種。對於整幅影象，統計出每一種(f1,f2)值出現的次數，然後排列成一個方陣，再用(f1,f2)出現的總次數將他們歸一化為出現的概率P(f1,f2)，由此產生的矩陣為灰度共生矩陣。θ方向上的間隔為d的灰度共生矩陣實際上是θ方向間隔為d的灰度變化量的聯合概率分佈。
2. 共生矩陣的計算
公式(1)中d表示畫素間隔，(k,l), (m,n)分別為原畫素和偏移後的畫素座標，其中k,m為縱座標，D為影象範圍[Image Processing, Analysis, and Machine Vision (Sonka 3rd Edition2007)]
舉例說明，假設原影象如圖1.a所示
對1.b中藍色字表示原畫素灰度值，紅字為偏移後像素灰度值。則對矩陣元素P0°,1 (0,0)表示1.a中在0°方向上（包括正和負方向）相距為1的(0,0)點對有兩對，考慮正負方向的加倍效果，P0°,1 (0,0)=4。同樣由於公式(1)對距離d定義的雙向性，使得灰度共生矩陣為對稱矩陣。
為了減小計算量，可將d定義為沿θ正方向。則(1)式變為
由1.a得到的新的灰度共生矩陣為
3. 共生矩陣計算紋理特徵
能量(Energy)：是灰度共生矩陣各元素值的平方和，是對影象紋理的灰度變化穩定程度的度量，反應了影象灰度分佈均勻程度和紋理粗細度。能量值大表明當前紋理是一種規則變化較為穩定的紋理。
熵(Entropy)：是影象包含資訊量的隨機性度量。當共生矩陣中所有值均相等或者畫素值表現出最大的隨機性時，熵最大；因此熵值表明了影象灰度分佈的複雜程度，熵值越大，影象越複雜。
最大概率(Maximum probability)：表示影象中出現次數最多的紋理特徵。
對比度(Contrast)：度量矩陣的值是如何分佈和影象中區域性變化的多少，反應了影象的清晰度和紋理的溝紋深淺。紋理的溝紋越深，反差越大，效果清晰；反之，對比值小，則溝紋淺，效果模糊。對公式(6)，典型的有κ=2，λ=1。
倒數差分矩(Inverse difference moment)：反映影象紋理的同質性，度量影象紋理區域性變化的多少。其值大則說明影象紋理的不同區域間缺少變化，區域性非常均勻。
相關性(Correlation)：自相關反應了影象紋理的一致性。如果影象中有水平方向紋理，則水平方向共生矩陣Correlation值大於其餘方向共生矩陣Correlation的值。它度量空間灰度共生矩陣元素在行或列方向上的相似程度，因此，相關值大小反映了影象中區域性灰度相關性。當矩陣元素值均勻相等時，相關值就大;相反，如果矩陣像元值相差很大則相關值小。
其中μx, μy為均值，σx, σy為標準差，計算公式如下
4. 演算法實現
#pragma once
#include<iostream>
#include <cassert>
#include <vector>
#include <iterator>
#include <functional>
#include <algorithm>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
typedef vector<vector<int> > VecGLCM;
typedef struct _GLCMFeatures
{
_GLCMFeatures()
: energy(0.0)
, entropy(0.0)
, contrast(0.0)
, idMoment(0.0)
{
}
double energy; // 能量
double entropy; // 熵
double contrast; // 對比度
double idMoment; // 逆差分矩， inverse difference moment
} GLCMFeatures;
class GLCM
{
public:
GLCM();
~GLCM();
public:
// 列舉灰度共生矩陣的方向
enum
{
GLCM_HORIZATION = 0, // 水平
GLCM_VERTICAL = 1, // 垂直
GLCM_ANGLE45 = 2, // 45度角
GLCM_ANGLE135 = 3 // 135度角
};
public:
// 計算灰度共生矩陣
// void calGLCM(IplImage* inputImg, VecGLCM& vecGLCM, int angle);
void calGLCM(cv::Mat &inputImg, VecGLCM &vecGLCM, int angle);
// 計算特徵值
void getGLCMFeatures(VecGLCM& vecGLCM, GLCMFeatures& features);
public:
// 初始化灰度共生矩陣
void initGLCM(VecGLCM& vecGLCM, int size = 16);
// 設定灰度劃分等級，預設值為 16
void setGrayLevel(int grayLevel) { m_grayLevel = grayLevel; }
// 獲取灰度等級
int getGrayLevel() const { return m_grayLevel; }
private:
// 計算水平灰度共生矩陣
void getHorisonGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
// 計算垂直灰度共生矩陣
void getVertialGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
// 計算 45 度灰度共生矩陣
void getGLCM45(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
// 計算 135 度灰度共生矩陣
void getGLCM135(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
private:
int m_grayLevel; // 將灰度共生矩陣劃分為 grayLevel 個等級
};#include "GLCM.h"
GLCM::GLCM() : m_grayLevel(16)
{
}
GLCM::~GLCM()
{
}
//==============================================================================
// 函式名稱: initGLCM
// 引數說明: vecGLCM,要進行初始化的共生矩陣,為二維方陣
// size, 二維矩陣的大小,必須與影象劃分的灰度等級相等
// 函式功能: 初始化二維矩陣
//==============================================================================
void GLCM::initGLCM(VecGLCM& vecGLCM, int size)
{
assert(size == m_grayLevel);
vecGLCM.resize(size);
for (int i = 0; i < size; ++i)
{
vecGLCM[i].resize(size);
}
for (int i = 0; i < size; ++i)
{
for (int j = 0; j < size; ++j)
{
vecGLCM[i][j] = 0;
}
}
}
//==============================================================================
// 函式名稱: getHorisonGLCM
// 引數說明: src,要進行處理的矩陣,源資料
// dst,輸出矩陣,計算後的矩陣，即要求的灰度共生矩陣
// imgWidth, 影象寬度
// imgHeight, 影象高度
// 函式功能: 計算水平方向的灰度共生矩陣
//==============================================================================
void GLCM::getHorisonGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight)
{
int height = imgHeight;
int width = imgWidth;
for (int i = 0; i < height; ++i)
{
for (int j = 0; j < width - 1; ++j)
{
int rows = src[i][j];
int cols = src[i][j + 1];
dst[rows][cols]++;
}
}
}
//==============================================================================
// 函式名稱: getVertialGLCM
// 引數說明: src,要進行處理的矩陣,源資料
// dst,輸出矩陣,計算後的矩陣，即要求的灰度共生矩陣
// imgWidth, 影象寬度
// imgHeight, 影象高度
// 函式功能: 計算垂直方向的灰度共生矩陣
//==============================================================================
void GLCM::getVertialGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight)
{
int height = imgHeight;
int width = imgWidth;
for (int i = 0; i < height - 1; ++i)
{
for (int j = 0; j < width; ++j)
{
int rows = src[i][j];
int cols = src[i + 1][j];
dst[rows][cols]++;
}
}
}
//==============================================================================
// 函式名稱: getGLCM45
// 引數說明: src,要進行處理的矩陣,源資料
// dst,輸出矩陣,計算後的矩陣，即要求的灰度共生矩陣
// imgWidth, 影象寬度
// imgHeight, 影象高度
// 函式功能: 計算45度的灰度共生矩陣
//==============================================================================
void GLCM::getGLCM45(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight)
{
int height = imgHeight;
int width = imgWidth;
for (int i = 0; i < height - 1; ++i)
{
for (int j = 0; j < width - 1; ++j)
{
int rows = src[i][j];
int cols = src[i + 1][j + 1];
dst[rows][cols]++;
}
}
}
//==============================================================================
// 函式名稱: getGLCM135
// 引數說明: src,要進行處理的矩陣,源資料
// dst,輸出矩陣,計算後的矩陣，即要求的灰度共生矩陣
// imgWidth, 影象寬度
// imgHeight, 影象高度
// 函式功能: 計算 135 度的灰度共生矩陣
//==============================================================================
void GLCM::getGLCM135(VecGLCM& src, VecGLCM& dst, int imgWidth, int imgHeight)
{
int height = imgHeight;
int width = imgWidth;
for (int i = 0; i < height - 1; ++i)
{
for (int j = 1; j < width; ++j)
{
int rows = src[i][j];
int cols = src[i + 1][j - 1];
dst[rows][cols]++;
}
}
}
//==============================================================================
// 函式名稱: calGLCM
// 引數說明: inputImg,要進行紋理特徵計算的影象,為灰度影象
// vecGLCM, 輸出矩陣,根據灰度影象計算出的灰度共生陣
// angle,灰度共生矩陣的方向,有水平、垂直、45度、135度四個方向
// 函式功能: 計算灰度共生矩陣
//==============================================================================
// void GLCM::calGLCM(IplImage* inputImg, VecGLCM& vecGLCM, int angle)
// {
// assert(inputImg->nChannels == 1);
// IplImage* src = NULL;
// src = cvCreateImage(cvGetSize(inputImg), IPL_DEPTH_32S, inputImg->nChannels);
// cvConvert(inputImg, src);
//
// int height = src->height;
// int width = src->width;
// int maxGrayLevel = 0;
// // 尋找最大畫素灰度最大值
// for (int i = 0; i < height; ++i)
// {
// for (int j = 0; j < width; ++j)
// {
// int grayVal = cvGetReal2D(src, i, j);
// if (grayVal > maxGrayLevel)
// {
// maxGrayLevel = grayVal;
// }
//
// }
// }// end for i
//
// ++maxGrayLevel;
// VecGLCM tempVec;
// // 初始化動態陣列
// tempVec.resize(height);
// for (int i = 0; i < height; ++i)
// {
// tempVec[i].resize(width);
// }
//
// if (maxGrayLevel > 16)//若灰度級數大於16，則將影象的灰度級縮小至16級，減小灰度共生矩陣的大小。
// {
// for (int i = 0; i < height; ++i)
// {
// for (int j = 0; j < width; ++j)
// {
// int tmpVal = cvGetReal2D(src, i, j);
// tmpVal /= m_grayLevel;
// tempVec[i][j] = tmpVal;
// }
// }
//
// if (angle == GLCM_HORIZATION) // 水平方向
// getHorisonGLCM(tempVec, vecGLCM, width, height);
// if (angle == GLCM_VERTICAL) // 垂直方向
// getVertialGLCM(tempVec, vecGLCM, width, height);
// if (angle == GLCM_ANGLE45) // 45 度灰度共生陣
// getGLCM45(tempVec, vecGLCM, width, height);
// if (angle == GLCM_ANGLE135) // 135 度灰度共生陣
// getGLCM135(tempVec, vecGLCM, width, height);
// }
// else//若灰度級數小於16，則生成相應的灰度共生矩陣
// {
// for (int i = 0; i < height; ++i)
// {
// for (int j = 1; j < width; ++j)
// {
// int tmpVal = cvGetReal2D(src, i, j);
// tempVec[i][j] = tmpVal;
// }
// }
//
// if (angle == GLCM_HORIZATION) // 水平方向
// getHorisonGLCM(tempVec, vecGLCM, width, height);
// if (angle == GLCM_VERTICAL) // 垂直方向
// getVertialGLCM(tempVec, vecGLCM, width, height);
// if (angle == GLCM_ANGLE45) // 45 度灰度共生陣
// getGLCM45(tempVec, vecGLCM, width, height);
// if (angle == GLCM_ANGLE135) // 135 度灰度共生陣
// getGLCM135(tempVec, vecGLCM, width, height);
// }
//
// cvReleaseImage(&src);
// }
void GLCM::calGLCM(cv::Mat & inputImg, VecGLCM & vecGLCM, int angle)
{
assert(inputImg.channels() == 1);
int height = inputImg.rows;
int width = inputImg.cols;
int maxGrayLevel = 0;
// 尋找最大畫素灰度最大值
for (int i = 0 ; i < height ; ++i)
{
for (int j = 0 ; j < width ; ++j)
{
int grayVal = inputImg.at<uchar>(i, j);
if (grayVal > maxGrayLevel)
{
maxGrayLevel = grayVal;
}
}
}
++maxGrayLevel;
VecGLCM tempVec;
tempVec.resize(height);
for (int i = 0 ; i < height ; ++i)
{
tempVec[i].resize(width);
}
if (maxGrayLevel > 16)
{
for (int i = 0 ; i < height ; ++i)
{
for (int j = 0 ; j < width ; ++j)
{
int tmpVal = inputImg.at<uchar>(i, j);
tmpVal /= m_grayLevel;
tempVec[i][j] = tmpVal;
}
}
if (angle == GLCM_HORIZATION)
{
getHorisonGLCM(tempVec, vecGLCM, width, height);
}
else if (angle == GLCM_VERTICAL)
{
getVertialGLCM(tempVec, vecGLCM, width, height);
}
else if (angle == GLCM_ANGLE45)
{
getGLCM45(tempVec, vecGLCM, width, height);
}
else if (angle == GLCM_ANGLE135)
{
getGLCM135(tempVec, vecGLCM, width, height);
}
}
else
{
for (int i = 0; i < height; ++i)
{
for (int j = 0; j < width; ++j)
{
int tmpVal = inputImg.at<uchar>(i, j);
tempVec[i][j] = tmpVal;
}
}
if (angle == GLCM_HORIZATION)
{
getHorisonGLCM(tempVec, vecGLCM, width, height);
}
else if (angle == GLCM_VERTICAL)
{
getVertialGLCM(tempVec, vecGLCM, width, height);
}
else if (angle == GLCM_ANGLE45)
{
getGLCM45(tempVec, vecGLCM, width, height);
}
else if (angle == GLCM_ANGLE135)
{
getGLCM135(tempVec, vecGLCM, width, height);
}
}
}
//==============================================================================
// 函式名稱: getGLCMFeatures
// 引數說明: vecGLCM, 輸入矩陣,灰度共生陣
// features,灰度共生矩陣計算的特徵值,主要包含了能量、熵、對比度、逆差分矩
// 函式功能: 根據灰度共生矩陣計算的特徵值
//==============================================================================
void GLCM::getGLCMFeatures(VecGLCM& vecGLCM, GLCMFeatures& features)
{
int total = 0;
for (int i = 0; i < m_grayLevel; ++i)
{
for (int j = 0; j < m_grayLevel; ++j)
{
total += vecGLCM[i][j]; // 求所有影象的灰度值的和
}
}
vector<vector<double> > temp;
temp.resize(m_grayLevel);
for (int i = 0; i < m_grayLevel; ++i)
{
temp[i].resize(m_grayLevel);
}
// 歸一化
for (int i = 0; i < m_grayLevel; ++i)
{
for (int j = 0; j < m_grayLevel; ++j)
{
temp[i][j] = (double)vecGLCM[i][j] / (double)total;
}
}
for (int i = 0; i < m_grayLevel; ++i)
{
for (int j = 0; j < m_grayLevel; ++j)
{
features.energy += temp[i][j] * temp[i][j];
if (temp[i][j] > 0)
features.entropy -= temp[i][j] * log(temp[i][j]); //熵
features.contrast += (double)(i - j)*(double)(i - j)*temp[i][j]; //對比度
features.idMoment += temp[i][j] / (1 + (double)(i - j)*(double)(i - j));//逆差矩
}
}
}
#include "GLCM.h"
int main()
{
IplImage* img = cvLoadImage("1.jpg", 0);
// cv::Mat src = cv::imread("1.jpg", 0);
GLCM glcm;
VecGLCM vec;
GLCMFeatures features;
glcm.initGLCM(vec);
// 水平
glcm.calGLCM(img, vec, GLCM::GLCM_HORIZATION);
glcm.getGLCMFeatures(vec, features);
// 垂直
glcm.calGLCM(img, vec, GLCM::GLCM_VERTICAL);
glcm.getGLCMFeatures(vec, features);
// 45 度
glcm.calGLCM(img, vec, GLCM::GLCM_ANGLE45);
glcm.getGLCMFeatures(vec, features);
// 135 度
glcm.calGLCM(img, vec, GLCM::GLCM_ANGLE135);
glcm.getGLCMFeatures(vec, features);
cout << "asm = " << features.energy << endl;
cout << "eng = " << features.entropy << endl;
cout << "Con = " << features.contrast << endl;
cout << "Idm = " << features.idMoment << endl;
system("pause");
return 0;
}