GLCM
#pragma once
#include <iostream>
//GLCM纹理特征
template<typename T>
class GLCM
{
public:
//构造函数
//参数：原始图像（二维数组），行数，列数，需要分级的数量，图像最小值，图像最大值
GLCM(T** const image, const int& rowNum, const int& colNum, const int& levelNum, const T& minValue, const T& maxValue):
m_colNum(colNum),m_rowNum(rowNum),m_levelNum(levelNum)
{
buildLevelImage(image, rowNum, colNum, levelNum, minValue, maxValue);
}
//析构函数
~GLCM() {
//释放比特化数组
for (int i = 0; i < m_rowNum; i++)
delete[] m_levelImage[i];
delete m_levelImage;
//释放GLCM矩阵
for (int i = 0; i < m_levelNum; i++)
delete[] m_GLCMMatrix[i];
delete m_GLCMMatrix;
//释放归一化GLCM矩阵
for (int i = 0; i < m_levelNum; i++)
delete[] m_normalizeGLCM[i];
delete m_normalizeGLCM;
}
public:
//建立比特化数组
void buildLevelImage(T** const image, const int& rowNum, const int& colNum, const int& levelNum, const T& minValue, const T& maxValue);
//计算GLCM矩阵&归一化
//输入步长，是否双向
void buildGLCM(const int stepI, const int stepJ, const bool bidirectional);
//计算特征
//对比度
float calCon();
//能量
float calEng();
//同质性
float calHom();
public:
//数组尺寸
const int m_rowNum;
const int m_colNum;
//分级数
const int m_levelNum;
//有效像元对
int m_valueCount;
//比特化数组
int** m_levelImage;
//GLCM数组
int** m_GLCMMatrix;
//归一化GLCM数组
float** m_normalizeGLCM;
};
template<typename T>
void GLCM<T>::buildLevelImage(T** const image, const int& rowNum, const int& colNum, const int& levelNum, const T& minValue, const T& maxValue)
{
//建立比特化数组
m_levelImage = new int* [m_rowNum];
for (int i = 0; i < m_rowNum; i++) {
m_levelImage[i] = new int[m_colNum];
for (int j = 0; j < m_colNum; j++) {
// 处理nodata值,用-1表示
if (image[i][j]<minValue || image[i][j]>maxValue)
m_levelImage[i][j] = -1;
else {
// 归一化
// 特殊情况：值刚好是上界，为了防止级数溢出，转成级数-1
if (image[i][j] == maxValue)
m_levelImage[i][j] = m_levelNum - 1;
else
m_levelImage[i][j] = static_cast<int>(((image[i][j] - minValue) / (maxValue - minValue)) * m_levelNum);
}
}
}
//测试代码，打印输出
//for (int i = 0; i < m_rowNum; i++) {
// for (int j = 0; j < m_colNum; j++) {
// std::cout << m_levelImage[i][j] << " ";
// }
// std::cout << std::endl;
//}
}
template<typename T>
void GLCM<T>::buildGLCM(const int stepI, const int stepJ, const bool bidirectional)
{
//初始化GLCM矩阵(赋值0)
m_GLCMMatrix = new int* [m_levelNum];
for (int i = 0; i < m_levelNum; i++) {
m_GLCMMatrix[i] = new int[m_levelNum];
for (int j = 0; j < m_levelNum; j++) {
m_GLCMMatrix[i][j] = 0;
}
}
//有效像元对数量
m_valueCount = 0;
//遍历图像，计算
for (int rowIndex = 0; rowIndex < m_rowNum; rowIndex++) {
for (int colIndex = 0; colIndex < m_colNum; colIndex++) {
//判断是否有效数据
if (m_levelImage[rowIndex][colIndex] == -1)
continue;
//对应位置是否有效（位置有效&值有效）
int tarPointI = rowIndex + stepI;
int tarPointJ = colIndex + stepJ;
if (tarPointI < 0 || tarPointI >= m_rowNum || tarPointJ < 0 || tarPointJ >= m_colNum)
continue;
if (m_levelImage[tarPointI][tarPointJ] == -1)
continue;
m_valueCount++;
int x = m_levelImage[rowIndex][colIndex];
int y = m_levelImage[tarPointI][tarPointJ];
m_GLCMMatrix[x][y] += 1;
}
}
//如果双向，再算一遍，只不过是减步长
for (int rowIndex = 0; rowIndex < m_rowNum; rowIndex++) {
for (int colIndex = 0; colIndex < m_colNum; colIndex++) {
//判断是否有效数据
if (m_levelImage[rowIndex][colIndex] == -1)
continue;
//对应位置是否有效（位置有效&值有效）
int tarPointI = rowIndex - stepI;
int tarPointJ = colIndex - stepJ;
if (tarPointI < 0 || tarPointI >= m_rowNum || tarPointJ < 0 || tarPointJ >= m_colNum)
continue;
if (m_levelImage[tarPointI][tarPointJ] == -1)
continue;
m_valueCount++;
int x = m_levelImage[rowIndex][colIndex];
int y = m_levelImage[tarPointI][tarPointJ];
m_GLCMMatrix[x][y] += 1;
}
}
//测试代码，打印GLCM
/*for (int i = 0; i < m_levelNum; i++) {
for (int j = 0; j < m_levelNum; j++)
std::cout << m_GLCMMatrix[i][j] << " ";
std::cout << std::endl;
}*/
//归一化
m_normalizeGLCM = new float* [m_levelNum];
for (int i = 0; i < m_levelNum; i++) {
m_normalizeGLCM[i] = new float[m_levelNum];
for (int j = 0; j < m_levelNum; j++)
m_normalizeGLCM[i][j] = m_GLCMMatrix[i][j] / static_cast<float>(m_valueCount);
}
//测试代码，打印归一化GLCM
//for (int i = 0; i < m_levelNum; i++) {
// for (int j = 0; j < m_levelNum; j++)
// std::cout << m_normalizeGLCM[i][j] << " ";
// std::cout << std::endl;
//}
}
template <typename T>
float GLCM<T>::calCon() {
float contrast = 0;
for (int i = 0; i < m_levelNum; i++) {
for (int j = 0; j < m_levelNum; j++)
contrast += (i - j) * (i - j) * m_normalizeGLCM[i][j];
}
return contrast;
}
template <typename T>
float GLCM<T>::calEng() {
float energy = 0;
for (int i = 0; i < m_levelNum; i++) {
for (int j = 0; j < m_levelNum; j++)
energy += m_normalizeGLCM[i][j] * m_normalizeGLCM[i][j];
}
return energy;
}
template <typename T>
float GLCM<T>::calHom() {
float homogeneity = 0;
for (int i = 0; i < m_levelNum; i++) {
for (int j = 0; j < m_levelNum; j++)
homogeneity += m_normalizeGLCM[i][j] / (1 + (i - j) * (i - j));
}
return homogeneity;
}
Moran‘s I
#pragma once
#include <iostream>
//莫兰指数
template<typename T>
class moran
{
public:
//构造函数
//参数：image原始图像，rowNum行数，columnNum列数，nodataValue无数据值
moran(T** const image, int rowNum, int columnNum, T nodataValue):
m_image(image),m_rowNum(rowNum),m_columnNum(columnNum),m_nodataValue(nodataValue)
{
cal_index();
cal_WeightMatrix();
}
//析构函数
~moran() {
//释放像元-权重矩阵索引
for (int rowIndex = 0; rowIndex < m_rowNum; rowIndex++)
delete[] m_matrixIndex[rowIndex];
delete[] m_matrixIndex;
//释放权重矩阵
for (int rowIndex = 0; rowIndex < m_pixelCount; rowIndex++)
delete[] m_weightMatrix[rowIndex];
delete[] m_weightMatrix;
//释放像元值-权重矩阵索引
delete[] m_valueIndex;
};
//计算栅格对应空间权重矩阵序号（映射矩阵）
void cal_index();
//计算空间权重矩阵（周围8像素相邻为1，否则为0）、
void cal_WeightMatrix();
//计算莫兰指数
double cal_moransI();
public:
//图像地址
T** m_image;
//图像行数
int m_rowNum;
//图像列数
int m_columnNum;
//图像无效值
T m_nodataValue;
//有效像元计数
int m_pixelCount;
//图像平均值
double m_imageMean;
//像元-权重矩阵索引
int** m_matrixIndex;
//像元值-权重矩阵索引
T* m_valueIndex;
//权重矩阵
int** m_weightMatrix;
//权重矩阵和
int m_weightMatrixSum;
};
template<typename T>
void moran<T>::cal_index()
{
//初始化有效像元计数
m_pixelCount = 0;
//初始化像元-权重矩阵索引（映射矩阵）
m_matrixIndex = new int* [m_rowNum];
for (int rowIndex = 0; rowIndex < m_rowNum; rowIndex++)
m_matrixIndex[rowIndex] = new int[m_columnNum];
//初始化像元值-权重矩阵索引
m_valueIndex = new T[m_columnNum * m_rowNum];
//计算映射&索引
for (int rowIndex = 0; rowIndex < m_rowNum; rowIndex++) {
for (int columnIndex = 0; columnIndex < m_columnNum; columnIndex++) {
if (m_image[rowIndex][columnIndex] == m_nodataValue)
continue;
else {
m_matrixIndex[rowIndex][columnIndex] = m_pixelCount;
m_valueIndex[m_pixelCount] = m_image[rowIndex][columnIndex];
m_pixelCount++;
}
}
}
}
template<typename T>
void moran<T>::cal_WeightMatrix()
{
//初始化权重矩阵
m_weightMatrix = new int* [m_pixelCount];
for (int rowIndex = 0; rowIndex < m_pixelCount; rowIndex++)
m_weightMatrix[rowIndex] = new int[m_pixelCount];
//赋初值
for (int rowIndex = 0; rowIndex < m_pixelCount; rowIndex++) {
for (int columnIndex = 0; columnIndex < m_pixelCount; columnIndex++)
m_weightMatrix[rowIndex][columnIndex] = 0;
}
//初始化权重矩阵和
m_weightMatrixSum = 0;
//初始化图像平均值
m_imageMean = 0;
//计算权重矩阵
for (int rowIndex = 0; rowIndex < m_rowNum; rowIndex++) {
for (int columnIndex = 0; columnIndex < m_columnNum; columnIndex++) {
//检查自身是否为空值
if (m_image[rowIndex][columnIndex] == m_nodataValue)
continue;
//更新平均值
m_imageMean += m_image[rowIndex][columnIndex];
//8相邻像元
for (int i = -1; i <= 1; i++) {
for (int j = -1; j <= 1; j++) {
//排除自身
if (i == 0 && j == 0)
continue;
int neighborRow = rowIndex + i;
int neighborColumn = columnIndex + j;
//检查是否范围内
if (neighborRow < 0 || neighborRow >= m_rowNum || neighborColumn < 0 || neighborColumn >= m_columnNum)
continue;
//检查是否为有效像元
if (m_image[neighborRow][neighborColumn] != m_nodataValue) {
//更新权重矩阵和
m_weightMatrixSum++;
//写入权重矩阵
int weightMatrixI = m_matrixIndex[rowIndex][columnIndex];
int weightMatrixJ= m_matrixIndex[neighborRow][neighborColumn];
m_weightMatrix[weightMatrixI][weightMatrixJ] = 1;
}
}
}
}
}
//计算平均值
m_imageMean /= m_pixelCount;
}
template<typename T>
double moran<T>::cal_moransI() {
//测试代码
//公式右上角累加
double rightup = 0;
for (int i = 0; i < m_pixelCount; i++) {
for (int j = 0; j < m_pixelCount; j++) {
rightup += m_weightMatrix[i][j] * (m_valueIndex[i] - m_imageMean) * (m_valueIndex[j] - m_imageMean);
}
}
//公式右下角累加
double rightdown = 0;
for (int i = 0; i < m_pixelCount; i++)
rightdown += (m_valueIndex[i] - m_imageMean) * (m_valueIndex[i] - m_imageMean);
//计算莫兰值
double moransI = 0;
moransI = m_pixelCount * rightup / (m_weightMatrixSum * rightdown);
return moransI;
}
Main函数测试
#include <iostream>
#include "GLCM.h"
#include "moran.h"
int main() {
const int numRows = 5;
const int numCols = 5;
const double minValue = 0.0;
const double maxValue = 1.0;
const int levelNum = 8;
//double testArray1[5][5] = { {0,0,0,0,0},{0,1,0,1,0}, {0,1,1,0,0}, {0,0,0,0,0}, {0,0,0,0,0} };
double** testArray = new double*[numRows];
srand((int)time(0));
//随机输入0-10
for (int i = 0; i < numRows; i++) {
testArray[i] = new double[numCols];
for (int j = 0; j < numCols; j++) {
// 生成随机值在minValue和maxValue之间
testArray[i][j] = minValue + (maxValue - minValue) * static_cast<double>(rand()) / RAND_MAX;
//testArray[i][j] = j>4?1:0;
//testArray[i][j] = i % 2;
//testArray[i][j] = testArray1[i][j];
}
}
//设置俩个空值试试
//testArray[0][0] = -999;
//estArray[3][4] = -999;
// 打印生成的二维数组
for (int i = 0; i < numRows; i++) {
for (int j = 0; j < numCols; j++) {
std::cout << testArray[i][j] << " ";
}
std::cout << std::endl;
}
//调用GLCM
GLCM<double> test_GLCM(testArray, numRows, numCols, levelNum, minValue, maxValue);
//计算步长0，1的GLCM
test_GLCM.buildGLCM(0, 1, false);
//计算同质度
float homogeneity = test_GLCM.calHom();
std::cout << "同质度：" << homogeneity << std::endl;
//计算对比度
float contrast = test_GLCM.calCon();
std::cout << "对比度：" << contrast << std::endl;
//计算能量
float Energy = test_GLCM.calEng();
std::cout << "能量：" << Energy << std::endl;
//计算莫兰指数
moran<double> test_moran(testArray,numRows,numCols,-9999);
double moransI = test_moran.cal_moransI();
std::cout << "莫兰指数：" << moransI << std::endl;
}
Comments NOTHING