基于视频模糊质量的评价（四）

2.5.2 基于梯度结构相似度的无参考模糊图像质量评价

这是结构相似度的一个扩展，详细的可以下载论文来进行查看。其算法的步骤如下：

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;
}

该算法的复杂度还是比较低的，且运算的结果也还行。