无参考图像的清晰度评价方法

在无参考图像的质量评价中，图像的清晰度是衡量图像质量优劣的重要指标，它能够较好的与人的主观感受相对应，图像的清晰度不高表现出图像的模糊。本文针对无参考图像质量评价应用，对目前几种较为常用的、具有代表性清晰度算法进行讨论分析，为实际应用中选择清晰度算法提供依据。
 

（1）Brenner 梯度函数

Brenner梯度函数是最简单的梯度评价函数，它只是简单的计算相邻两个像素灰度差的平方，该函数定义如下：

        其中：f(x,y) 表示图像f对应像素点(x,y)的灰度值，D(f)为图像清晰度计算结果(下同)。

 

（2）Tenengrad 梯度函数

Tenengrad 梯度函数采用Sobel算子分别提取水平和垂直方向的梯度值，基与Tenengrad 梯度函数的图像清晰度定义如下：

       G(x,y) 的形式如下： 

  其中：T是给定的边缘检测阈值，Gx和Gy分别是像素点(x,y)处Sobel水平和垂直方向边缘检测算子的卷积，建议使用以下的Sobel算子模板来检测边缘：

 

（3）Laplacian 梯度函数

Laplacian 梯度函数与Tenengrad梯度函数基本一致，用Laplacian算子替代Sobel算子即可，该算子定义如下：

        因此基于Laplacian 梯度函数的图像星清晰度的定义如下：

        其中G(x,y)是像素点(x,y)处Laplacian算子的卷积。

 

（4）SMD（灰度方差）函数

        当完全聚焦时，图像最清晰，图像中的高频分量也最多，故可将灰度变化作为聚焦评价的依据，灰度方差法的公式如下：

 

（5）SMD2 （灰度方差乘积）函数

        灰度差分评价函数具有较好的计算性能，但其缺点也很明显，即在焦点附近灵敏度不高，即该函数在极值点附近过于平坦，从而导致聚焦精度难以提高。在文章《一种快速高灵敏度聚焦评价函数》中提出了一种新的评价函数，称之为灰度方差乘积法，即对每一个像素领域两个灰度差相乘后再逐个像素累加，该函数定义如下：

 

（6）方差函数

        因为清晰聚焦的图像有着比模糊图像更大的灰度差异，可以将方差函数作为评价函数：

        

  其中：为整幅图像的平均灰度值，该函数对噪声比较敏感，图像画面越纯净，函数值越小。

 

（7）能量梯度函数

        能量梯度函数更适合实时评价图像清晰度，该函数定义如下：

 

（8）Vollath函数

        Vollath函数定义如下：

 其中：为整幅图像的平均灰度值，M和N分别为图像宽和高。
  

（9）熵函数

       基于统计特征的熵函数是衡量图像信息丰富程度的一个重要指标，有信息论可知，一幅图像 f 的信息量是由该图像的信息熵 D(f) 来度量：

 其中：Pi 是图像中灰度值为i的像素出现的概率，L为灰度级总数（通常取值256）。根据Shannon信息论，熵最大时信息量最多。将此原理应用到对焦过程，D(f)越大则图像越清晰。熵函数灵敏度不高，依据图像内容不同容易出现与真实情况相反的结果。

 

（10） EAV点锐度算法函数

        徐贵力、张霞等提出了一种基于边缘锐度的算法用于评价图像的清晰度。通过统计图像某一边缘方向的灰度变化情况来评价。计算公式如下：

 其中：df/dx为边缘法向的灰度变化率，f(b) - f(a)为该方向的总体灰度变化。该算法只对图像的特定边缘区域做统计，能否代表整幅图像的清晰度仍有疑问，此外计算前需人工选定边缘区域，不便实现程序运算的自动化，因为王鸿南等在论文 图像清晰度评价方法研究 中对上述算法进行了改进，改进如下：

        a)  将针对边缘的梯度计算改为逐个像素领域梯度的计算，以便算法能对图像的整体进行评价，并使算法实现自动化。

        b)  对方格像素 8 领域的灰度变化进行距离加权，水平和垂直方向的权重为1，而45度和135度方向的权重为 。

        c)  对计算结果按图像的大小进行规格化，以便于图像的对比。

经过以上三步改进后的点锐度算法为：

 其中：M和N为图像的行数和列数。

 

（11）Reblur 二次模糊
         如果一幅图像已经模糊了，那么再对它进行一次模糊处理，高频分量变化不大；但如果原图是清楚的，对它进行一次模糊处理，则高频分量变化会非常大。因此可以通过对待评测图像进行一次高斯模糊处理，得到该图像的退化图像，然后再比较原图像和退化图像相邻像素值的变化情况，根据变化的大小确定清晰度值的高低，计算结果越小表明图像越清晰，反之越模糊。这种思路可称作基于二次模糊的清晰度算法，其算法简化流程如下图：

 

 

 

（12）NRSS 梯度结构相似度
          Wang等利用人类视觉系统（HVS）非常适于提取目标的结构信息的特点，提出了图像结构相似度概念（SSIM），认为只要能计算目标结构信息的变化，就能够得到感知图像失真值。杨春玲等基于此思路，将该方法引入到计算全参考图像的清晰度评价中，认为图像的清晰度可以使用目标图像与参考图像间的结构相似度来表示，而图像间的结构相似度包含以下三个部分的比较：

 而C1、C2和C3 是为了避免分母为0而设的常数。图像的结构相似度由下式计算可得：

 为简单起见可以令 
        谢小甫等进一步改进了杨春玲等的方法，根据结构相似度的相关思想结合人烟视觉系统的相关特点，设计了无参考图像清晰度的评价指标（NRSS），计算方法如下：

 

 

   （a）为待评价图像构造参考图像。定义待评价图像为I，而参考图像 ，即对待评价图像I进行低通滤波得到参考图像 。实验表明，基于圆盘模型的均值滤波器和高斯模型的平滑滤波器都可以取得较好的效果，为了更好的与成像系统匹配，建议采用7x7大小且 的高斯平滑滤波器。在需要实时处理的工程应用中7x7均值滤波器并不会是评价效果下降很大。
   （b）提取图像 I 和 的梯度信息。利用人眼对水平和垂直方向的边缘信息最为敏感的特性，使用Sobel算子分别提取水平和垂直方向的边缘信息，定义 I 和 的梯度图像是G 和 。
   （c）找出梯度图像 G 中梯度信息最丰富的 N 个图像块。将图像G划分为 8x8 的小块，块间的步长为4，即相邻块有50%重叠，这是为了避免丢失重要的边缘。计算每块的方差，方差越大说明梯度信息越丰富，找出其中方差最大的N块，记为，对应的中的对应块定义为  。N的值大小直接影响评价结果，同时也影响算法运行时间。在后面的实验中设 N = 64。
   （d）计算结构清晰度NRSS。先计算每个xi和 yi的结构相似度SSIM(xi, yi)，其中SSIM计算方法参见前面的定义，则图像的无参考结构清晰度定义为：

 

（13）FFT 图像变换域

      待写！

 

（14）No-Reference Perceptual Quality Assessment of JPEG Compressed Images

        在这篇文章中，作者【Zhou Wang】等针对JPEG压缩图片提出了一种新的无参图像质量评价方法。
         JPEG图片是基于8x8块的DCT变换的编码技巧，它是有损的因为对DCT变换系数做量化的时候会产生量化误差。量化就会导致模糊和块效应。模糊主要是因为丢失了高频的DCT系数。块效应是由于块边界的不连续性，因为每个分块的量化是独立的。

        我们用 f(x, y) 表示一幅图片，图片尺寸为 MxN，计算跨越每个水平线的信号差：

        首先计算块效应，块效应的定义就是左右跨越边界的信号差的平均值：

        然后计算块内信号差的平均值：

       再计算zero-crossing（ZC）率，ZC是边界跨零的意思，也就是说相邻两个点的值的乘积为负数，也就是一正一负，因此对于[1, N - 2]范围内的y，定义如下变量：

 

 于是水平方向的ZC率定义如下：

         同理，我们可以计算垂直方向的几个指标值 。最后得到这几个指标的水平和垂直方向的平均值：

 

 有很多方式把这几个指标联系起来组成一个质量评价模型。此处我们采用如下图像质量定义：

 其中 是从大量实验中提炼出来的模型参数。本文中所采用的参数值如下：

 

 

（15）No-Reference Image Quality Assessment forJPEG/JPEG2000 Coding

         这篇文章的作者在前面那篇文章的基础上，重新定义了新的质量指标：

 其实 S 就是在（14）中已经得到的质量评价值。

 

（16）No-Reference Image Quality Assessment  using Blur and Noise
     图像质量受很多因素影响，例如：亮度、对比度、色调、边界、噪声、模糊等。在本文中，我们假定噪声和模糊是影响图像质量最重要的两个因素。简单起见，只对彩色图像的亮度分量做模糊和噪声监测。本文的图像质量评价算法框架图如下：

    A）模糊检测
      模糊估计分为两个步骤：首先是边缘检测，然后是模糊确定。此处模糊估计是通过计算当前像素点与领域内像素点均值之差来确定。我们用f(x,y) 表示图片，其中。定义水平绝对差如下：

 整个图片的水平绝对差的均值为:

 如果当前像素点的  则该像素点就是一个候选的边缘点. 如果比它水平方向两个相邻的点 都大，则该像素点就被确认为一个边缘点。边缘点的判断总结如下： 

 

  接下来我们检测边缘点是否模糊。定义：

同理，按照以上的步骤我们可以计算垂直方向的值 。两者之大者称作Inverse Blurriness，用于最终的模糊判定依据。

 低于阈值ThB的Inverse Blurriness 被认为是模糊的。实验测试表明此处的阈值ThB取值0.1。最后，边缘模糊的均值和比率为：

  

    B）噪点检测
      因为沿边缘的噪点视觉上不明显，因此我们只检测边缘之外的噪点。边缘检测会被噪点影响，因此在检测边缘之前做一个噪点滤波的预处理。在本文中，我们应用均值滤波来消除噪点。均值滤波之后的图像g(x,y)为：

候选的噪点估计如下：

同理可以在垂直方向计算对应的值。然后得到候选的噪点是：

其中N_cand(x,y)表示候选噪点，它在边缘区域为0。

噪点均值和比率为：

其中Sum_Noise是N(x,y)之和，Noise_cnt是噪点总数目。

   C）噪点和模糊的组合
      此处我们的图像质量评价指标定义如下：

其中w1、w2、w3、w4是权值。通过线性回归分析获取这些权值。本文中这些权值为：

 

 

部分源码

void getImageQualityAssessment(cv::Mat img) {

    //可以直接根据灰度图的方差 做判断

    //Tenengrad 梯度
    cv::Mat imageSobel;
    Sobel(img, imageSobel, CV_16U, 1, 1);
    double meanValue = 0.0;
    meanValue = mean(imageSobel)[0];

    //Laplacian
    cv::Mat dst;
    //Laplace(f) = \dfrac{\partial^{2} f}{\partial x^{2}} + \dfrac{\partial^{2} f}{\partial y^{2}}
    int kernel_size = 3;
    int ddepth = CV_16U;
    Laplacian(img, dst, ddepth, kernel_size, BORDER_DEFAULT);
    //cv::imshow("laplac", dst);
    Mat tmp_m, tmp_sd;
    double m = 0, sd = 0;
    meanStdDev(dst, tmp_m, tmp_sd);
    m = tmp_m.at<double>(0, 0);
    sd = tmp_sd.at<double>(0, 0);
    cout << "laplacian......... Mean: " << m << " ,laplacian StdDev: " << sd << endl;

}

 

#include "opencv2/opencv.hpp"
#include "iomanip"

using namespace cv;
using namespace std;

unsigned long Tenengrad_measure(Mat& image)
{
  unsigned long score = 0;
  for (int x = 0; x < image.rows; ++x)
  {
    uchar* ptr = image.ptr<uchar>(x);
#pragma omp parallel for
    for (int y = 0; y < image.cols*image.channels(); ++y)
    {
      score += ptr[y];
    }
  }
  return score;
}

double Tenengrad(Mat& roi)
{
  Mat smooth_image;
  blur(roi, smooth_image, Size(3, 3));
  Mat sobel_x_mat, sobel_y_mat;
  Sobel(smooth_image, sobel_x_mat, -1, 1, 0, 0);
  Sobel(smooth_image, sobel_y_mat, -1, 0, 1, 3);

  convertScaleAbs(sobel_x_mat, sobel_x_mat);
  convertScaleAbs(sobel_y_mat, sobel_y_mat);

  Mat pow_x;
  Mat pow_y;
  pow(sobel_x_mat, 2, pow_x);
  pow(sobel_y_mat, 2, pow_y);

  Mat weighted_mat;
  addWeighted(pow_x, 1, pow_y, 1, 0, weighted_mat);

  double score = Tenengrad_measure(weighted_mat);

  return score;
}

double Brenner(Mat& image)
{
  double score = 0.0;
  for (int i = 0; i < image.rows; ++i)
  {
    uchar *ptr = image.ptr<uchar>(i);
#pragma omp parallel for
    for (int j = 0; j < (image.cols - 2)*image.channels(); ++j)
    {
      score += (ptr[j + 2] - ptr[j])*(ptr[j + 2] - ptr[j]);
    }
  }
  return score;
}

double SMD2(Mat& image)
{
  unsigned long score = 0;
  for (int i = 0; i < image.rows - 1; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
    uchar* ptr_1 = image.ptr<uchar>(i + 1);
#pragma omp parallel for
    for (int j = 0; j < (image.cols - 1)*image.channels(); ++j)
    {
      score += (ptr[j] - ptr[j + 1]) * (ptr[j] - ptr_1[j]);
    }
  }
  return score;
}

double Energy(Mat& image)
{
  unsigned long score = 0;
  for (int i = 0; i < image.rows - 1; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
    uchar* ptr_1 = image.ptr<uchar>(i + 1);
#pragma omp parallel for
    for (int j = 0; j < (image.cols - 1) * image.channels(); ++j)
    {
      score += pow(ptr[j + 1] - ptr[j], 2) + pow(ptr_1[j] - ptr[j], 2);
    }
  }
  return score;
}

double Jpeg(Mat& image)
{
  //horizontal calculate
  auto b_h = 0.0;
  for (int i = 1; i < floor(image.rows / 8) - 1; ++i)
  {
    uchar* ptr = image.ptr<uchar>(8 * i);
    uchar* ptr_1 = image.ptr<uchar>(8 * i + 1);
#pragma omp parallel for
    for (int j = 1; j < image.cols; ++j)
    {
      b_h += abs(ptr_1[j] - ptr[j]);
    }
  }
  b_h *= 1 / (image.cols*(floor(image.rows / 8) - 1));

  auto a_h = 0.0;
  for (int i = 1; i < image.rows - 1; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
    uchar* ptr_1 = image.ptr<uchar>(i + 1);
#pragma omp parallel for
    for (int j = 1; j < image.cols; ++j)
    {
      a_h += abs(ptr_1[j] - ptr[j]);
    }
  }
  a_h = (a_h * 8.0 / (image.cols * (image.rows - 1)) - b_h) / 7;

  auto z_h = 0.0;
  for (int i = 1; i < image.rows - 2; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
    uchar* ptr_1 = image.ptr<uchar>(i + 1);
#pragma omp parallel for
    for (int j = 1; j < image.cols; ++j)
    {
      z_h += (ptr_1[j] - ptr[j]) * (ptr_1[j + 1] - ptr[j]) > 0 ? 0 : 1;
    }
  }
  z_h *= 1.0 / (image.cols* (image.rows - 2));

  //vertical calculate
  auto b_v = 0.0;
  for (int i = 1; i < image.rows; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
#pragma omp parallel for
    for (int j = 1; j < floor(image.cols / 8) - 1; ++j)
    {
      b_v += abs(ptr[8 * j + 1] - ptr[8 * j]);
    }
  }
  b_v *= 1.0 / (image.rows*(floor(image.cols / 8) - 1));

  auto a_v = 0.0;
  for (int i = 1; i < image.rows; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
#pragma omp parallel for
    for (int j = 1; j < image.cols - 1; ++j)
    {
      a_v += abs(ptr[j + 1] - ptr[j]);
    }
  }
  a_v = (a_v * 8.0 / (image.rows * (image.cols - 1)) - b_v) / 7;

  auto z_v = 0.0;
  for (int i = 1; i < image.rows; ++i)
  {
    uchar* ptr = image.ptr<uchar>(i);
#pragma omp parallel for
    for (int j = 1; j < image.cols - 2; ++j)
    {
      z_v += (ptr[j + 1] - ptr[j]) * (ptr[j + 2] - ptr[j + 1]) > 0 ? 0 : 1;
    }
  }
  z_v *= 1.0 / (image.rows* (image.cols - 2));

  
  auto B = (b_v + b_h) / 2;
  auto A = (a_h + a_v) / 2;
  auto Z = (z_h + z_v) / 2;
  auto alpha = -245.9, beta = 261.9, gamma1 = -0.024, gamma2 = 0.016, gamma3 = 0.0064;

  auto S = alpha + beta*pow(B, gamma1)*pow(A, gamma2)*pow(Z, gamma3);

  return S;
}

double Jpeg2(Mat& image)
{
  double s = Jpeg(image);
  double ss = 4.0 / (1.0 + exp(-1.0217*(s - 3))) + 1.0;
  return ss;
}

int main()
{
  VideoCapture capture;
  if (!capture.open("C://Users//sp_zh//Desktop//1.avi"))
  {
    return -1;
  }
  int frames = capture.get(CV_CAP_PROP_FRAME_COUNT);
  int index = 0;
  vector<double> focuss;
  while (frames - index > 0)
  {
    Mat image;
    capture >> image;

    Rect roi_rect(1221, 409, 493, 1267);
    Mat roi = image(roi_rect);
    //double score = Tenengrad(roi);
    //double score = Brenner(roi);
    //double score = SMD2(roi);
    //double score = Energy(roi);
    //double score = Jpeg(roi);
    double score = Jpeg2(roi);
    char text[16];
    sprintf(text, "%lf", score);
    putText(image, text, Point(100, 100), FONT_HERSHEY_COMPLEX, 1.5, Scalar(255, 0, 0));
    imshow("1", image);
    waitKey(20);
    index++;
    focuss.push_back(score);
  }
  return 0;
}