1.先将原始的图像使用11 x 11且方差为1的高斯滤波器进行滤波
2.在接着对模糊后的图像进行半径为4的菱形的Sobel算子的边缘提取
3.将提取后的内容进行膨胀
4.在膨胀图像的情况下判断膨胀模块，并在膨胀模块下计算其GSIM值（相关的膨胀子块可以查看论文，论文不难理解）
5.最后将所有的GSIM值加起来除于膨胀模块的个数，便得到模糊值。
相关代码的如下，有点凌乱，且没有加保护措施：
double DefGSIM(Mat frame) { double si1,sj1,si2,sj2,Gx,Gy,temp,sum=0,result,count =0; int i,j; int scale = 1; int delta = 0; int ddepth = CV_16S; double k1=0.01; int L=255; double c1=pow((k1*L), 2); Mat gray,Gaussian,sobel,dilate_result; cvtColor(frame, gray, CV_BGR2GRAY); IplImage *src = &(IplImage(gray)); uchar *src_data = (uchar*)src->imageData; GaussianBlur(gray, Gaussian, Size(11, 11), 1); IplImage *gua = &(IplImage(Gaussian)); uchar *gua_data = (uchar*)src->imageData; //Sobel operation /// 创建 grad_x 和 grad_y 矩阵 Mat grad_x, grad_y; Mat abs_grad_x, abs_grad_y; /// 求 X方向梯度 //Scharr( src_gray, grad_x, ddepth, 1, 0, scale, delta, BORDER_DEFAULT ); Sobel( Gaussian, grad_x, ddepth, 1, 0, 3, scale, delta, BORDER_DEFAULT ); convertScaleAbs( grad_x, abs_grad_x ); /// 求Y方向梯度 //Scharr( src_gray, grad_y, ddepth, 0, 1, scale, delta, BORDER_DEFAULT ); Sobel( Gaussian, grad_y, ddepth, 0, 1, 3, scale, delta, BORDER_DEFAULT ); convertScaleAbs( grad_y, abs_grad_y ); /// 合并梯度(近似) addWeighted( abs_grad_x, 0.5, abs_grad_y, 0.5, 0, sobel ); Mat element = getStructuringElement( MORPH_RECT, Size(4,4 ) ); /// 膨胀操作 dilate( sobel, dilate_result, element ); IplImage *img = &(IplImage(dilate_result)); int height = img->height; int width = img->width; int step = img->widthStep/sizeof(uchar); uchar *data = (uchar*)img->imageData; Mat gsim = Mat::zeros(height, width, CV_64F); for (i = 1; i<height-1; i++) { for(j = 1; j<width-1; j++) { if ((data[i*step+j] == 1) && (data[(i-1)*step+j-1] == 1) && (data[(i-1)*step+j] == 1) && (data[(i+1)*step+j] == 1) && (data[i*step+j-1] == 1) && (data[i*step+j+1] == 1) && (data[(i+1)*step+j-1] == 1) && (data[(i+1)*step+j] == 1) && (data[(i+1)*step+j+1] == 1)) { si1 = 3*src_data[(i-1)*step+j-1] + 10*src_data[i*step+j-1] + 3*src_data[(i+1)*step+j-1] - 3*src_data[(i-1)*step+j+1] - 10*src_data[i*step+j+1] - 3*src_data[(i+1)*step+j+1]; sj1 = 3*src_data[(i-1)*step+j-1] + 10*src_data[(i-1)*step+j] + 3*src_data[(i-1)*step+j+1] - 3*src_data[(i+1)*step+j-1] - 10*src_data[(i+1)*step+j] - 3*src_data[(i+1)*step+j+1]; Gx = abs(si1) + abs(sj1); si2 = 3*gua_data[(i-1)*step+j-1] + 10*gua_data[i*step+j-1] + 3*gua_data[(i+1)*step+j-1] - 3*gua_data[(i-1)*step+j+1] - 10*gua_data[i*step+j+1] - 3*gua_data[(i+1)*step+j+1]; sj2 = 3*gua_data[(i-1)*step+j-1] + 10*gua_data[(i-1)*step+j] + 3*gua_data[(i-1)*step+j+1] - 3*gua_data[(i+1)*step+j-1] - 10*gua_data[(i+1)*step+j] - 3*gua_data[(i+1)*step+j+1]; Gy = abs(si2) + abs(sj2); temp = (2*Gx*Gy + c1)/(Gx*Gx+Gy*Gy+c1); gsim.at<double>(i,j) = temp ; count = count + 1; } } } for (int i = 0; i < height;i++) { uchar* data = gsim.ptr<uchar>(i); for (int j = 0; j < width; j++) { sum = data[j] + sum; } } result = sum/count; return result; }该算法的复杂度还是比较低的，且运算的结果也还行。