使用灰度共生矩阵实现指纹分割 C++

#include <iostream> #include "highgui.h" #include "cv.h" #define GLCM_DIS 2 //灰度共生矩阵的统计距离 #define GLCM_CLASS 8 //计算灰度共生矩阵的图像灰度值等级化 #define GLCM_ANGLE_HORIZATION 0 //水平 #define GLCM_ANGLE_VERTICAL 1 //垂直 #define GLCM_ANGLE_DIGONAL_45 2 //45度对角 #define GLCM_ANGLE_DIGONAL_135 3 //135度对角 int threshold = 15; #pragma region calGLCM int GetMax(CvMat* bWavelet,int& max_value,int& min_value) { int i,j; int width,height; int max = 0; int min = 0; width = bWavelet->width; height = bWavelet->height; for (i = 0;i<height;i++ ) { for(j = 0;j<width;j++) { if (CV_MAT_ELEM(*bWavelet,float,i,j)> max) max = CV_MAT_ELEM(*bWavelet,float,i,j); if(CV_MAT_ELEM(*bWavelet,float,i,j) < min) min = CV_MAT_ELEM(*bWavelet,float,i,j); } } max_value = max; min_value = min; return 1; } int calGLCM(CvMat* bWavelet,int angleDirection,double* featureVector) { int i,j; int width,height; int min,max; 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; //GetMax(bWavelet,max,min); //printf("max = %d,min = %d ",max,min); //system("pause"); int relative_value = 0; for(i = 0;i < height;i++){ for(j = 0;j < width;j++){ histImage[i * width + j] = (int)(CV_MAT_ELEM(*bWavelet,float,i,j) * GLCM_CLASS / 256); //relative_value = CV_MAT_ELEM(*bWavelet,float,i,j) - min; //histImage[i * width + j] = (int)(relative_value * GLCM_CLASS / (max-min+1)); //printf("%d ",histImage[i * width + j]); } } //初始化共生矩阵 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]++; } } } } //135度对角方向 else if(angleDirection == GLCM_ANGLE_DIGONAL_135) { 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]++; } } } } //45度对角方向 else if(angleDirection == GLCM_ANGLE_DIGONAL_45) { 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; } #pragma endregion calcGLCM double calc_variance(double a[4]) { double sum; double mean; double variance; sum = a[0]+a[1]+a[2]+a[3]; mean = sum/4.0; variance = (a[0]-mean)*(a[0]-mean) + (a[1]-mean)*(a[1]-mean) + (a[2]-mean)*(a[2]-mean)+(a[3]-mean)*(a[3]-mean); variance = variance/4.0; return sqrt(variance); } void main() { clock_t start,finish; int WIDTH,HEIGHT; threshold = 8.5; IplImage* src_img = cvLoadImage("1.bmp",0); //load the picture WIDTH = src_img->width; HEIGHT = src_img->height; int step = src_img->widthStep; //the width of each row printf("width: %d step: %d depth: %d channels: %d",WIDTH,step,src_img->depth,src_img->nChannels); IplImage* msk_img = cvCreateImage(cvSize(WIDTH,HEIGHT),8,1); cvZero(msk_img); IplImage* dst_img = cvCreateImage(cvSize(WIDTH,HEIGHT),8,1); cvZero(dst_img); cvNot(dst_img,dst_img); start = clock(); //cvSmooth(src_img,src_img,CV_GAUSSIAN); //cvSmooth(src_img,src_img,CV_MEDIAN,2); float sum; cvNamedWindow("src"); cvShowImage("src",src_img); CvMat* cwindow = cvCreateMat(5,4,CV_32FC1); //store the temperate 3*3 matrix //CvMat* cwindow = cvCreateMat(5,5,CV_32FC1); //store the temperate 3*3 matrix double features[10]; double entropy[4] = {0}; double energy[4] = {0}; double contrast[4] = {0}; double homogenity[4] = {0}; double variance; uchar* gray_ptr = NULL; uchar* msk_ptr = NULL; gray_ptr = (uchar*) (src_img->imageData); msk_ptr = (uchar*) (msk_img->imageData); int x,y; for(y = 2; y<src_img->height-2; y = y+2) { for ( x = 2; x<src_img->width-2; x = x+2) { /*CV_MAT_ELEM(*cwindow,float,0,0) = gray_ptr[(y-2)*step+(x-2)]; CV_MAT_ELEM(*cwindow,float,0,1) = gray_ptr[(y-2)*step+x-1]; CV_MAT_ELEM(*cwindow,float,0,2) = gray_ptr[(y-2)*step+x]; CV_MAT_ELEM(*cwindow,float,0,3) = gray_ptr[(y-2)*step+(x+1)]; CV_MAT_ELEM(*cwindow,float,0,4) = gray_ptr[(y-2)*step+(x+2)]; CV_MAT_ELEM(*cwindow,float,1,0) = gray_ptr[(y-1)*step+(x-2)]; CV_MAT_ELEM(*cwindow,float,1,1) = gray_ptr[(y-1)*step+x-1]; CV_MAT_ELEM(*cwindow,float,1,2) = gray_ptr[(y-1)*step+ x]; CV_MAT_ELEM(*cwindow,float,1,3) = gray_ptr[(y-1)*step+(x+1)]; CV_MAT_ELEM(*cwindow,float,1,4) = gray_ptr[(y-1)*step+x+2]; CV_MAT_ELEM(*cwindow,float,2,0) = gray_ptr[y*step+(x-2)]; CV_MAT_ELEM(*cwindow,float,2,1) = gray_ptr[y*step+(x-1)]; CV_MAT_ELEM(*cwindow,float,2,2) = gray_ptr[y*step+x]; CV_MAT_ELEM(*cwindow,float,2,3) = gray_ptr[y*step+(x+1)]; CV_MAT_ELEM(*cwindow,float,2,4) = gray_ptr[y*step+(x+2)]; CV_MAT_ELEM(*cwindow,float,3,0) = gray_ptr[(y+1)*step+(x-2)]; CV_MAT_ELEM(*cwindow,float,3,1) = gray_ptr[(y+1)*step+x-1]; CV_MAT_ELEM(*cwindow,float,3,2) = gray_ptr[(y+1)*step+x]; CV_MAT_ELEM(*cwindow,float,3,3) = gray_ptr[(y+1)*step+(x+1)]; CV_MAT_ELEM(*cwindow,float,3,4) = gray_ptr[(y+1)*step+(x+2)]; CV_MAT_ELEM(*cwindow,float,4,0) = gray_ptr[(y+2)*step+(x-2)]; CV_MAT_ELEM(*cwindow,float,4,1) = gray_ptr[(y+2)*step+x-1]; CV_MAT_ELEM(*c