有没有办法检测图像是否模糊？

Question

Closed. This question needs to be more focused. It is not currently accepting answers.

---

想要改进这个问题吗？通过编辑这篇文章来更新问题，使其仅关注一个问题。

3 年前已关闭。

Answer 1

估计图像清晰度的另一种非常简单的方法是使用拉普拉斯（或 LoG）滤波器并简单地选择最大值。如果您预计会有噪声（即选择第 N 个最高对比度而不是最高对比度），那么使用像 99.9% 分位数这样的稳健测量可能会更好。如果您预计图像亮度会发生变化，则还应该包括一个预处理步骤来标准化图像亮度/对比度（例如直方图均衡）。

我已经在 Mathematica 中实现了 Simon 的建议和这个建议，并在一些测试图像上进行了尝试：

第一个测试使用具有不同内核大小的高斯滤波器模糊测试图像，然后计算模糊图像的 FFT 并取 90% 最高频率的平均值：

testFft[img_] := Table[
  (
   blurred = GaussianFilter[img, r];
   fft = Fourier[ImageData[blurred]];
   {w, h} = Dimensions[fft];
   windowSize = Round[w/2.1];
   Mean[Flatten[(Abs[
       fft[[w/2 - windowSize ;; w/2 + windowSize, 
         h/2 - windowSize ;; h/2 + windowSize]]])]]
   ), {r, 0, 10, 0.5}]

得出对数图：

5行代表5张测试图像，X轴代表高斯滤波器半径。图形逐渐减小，因此 FFT 是衡量锐度的一个很好的方法。

这是“最高 LoG”模糊估计器的代码：它只是应用 LoG 滤波器并返回滤波器结果中最亮的像素：

testLaplacian[img_] := Table[
  (
   blurred = GaussianFilter[img, r];
   Max[Flatten[ImageData[LaplacianGaussianFilter[blurred, 1]]]];
   ), {r, 0, 10, 0.5}]

对数图的结果：

这里未模糊图像的分布要好一些（2.5 vs 3.3），主要是因为该方法仅使用图像中最强的对比度，而 FFT 本质上是整个图像的平均值。函数也减少得更快，因此设置“模糊”阈值可能更容易。

Answer 2

是的。计算快速傅立叶变换并分析结果。傅里叶变换告诉您图像中存在哪些频率。如果高频数量较少，则图像会模糊。

定义术语“低”和“高”取决于您。

编辑：

如评论中所述，如果您想要一个代表给定图像的模糊度的单个浮点数，则必须计算出合适的指标。

nikie 的回答提供了这样一个指标。将图像与拉普拉斯核进行卷积：

   1
1 -4  1
   1

并在输出上使用稳健的最大度量来获取可用于阈值化的数字。在计算拉普拉斯算子之前，尽量避免对图像进行过度平滑，因为您只会发现平滑后的图像确实是模糊的:-)。

Answer 3

在使用自动对焦镜头的一些工作中，我遇到了这套非常有用的算法 检测图像焦点。它是在 MATLAB 中实现的，但大多数功能很容易通过 filter2D。

它基本上是许多焦点测量算法的调查实现。如果您想阅读原始论文，代码中提供了算法作者的参考资料。 Pertuz 等人 2012 年发表的论文。 焦点分析焦点形状测量操作符 (SFF) 对所有这些测量及其性能（在应用于 SFF 的速度和准确性方面）进行了很好的回顾。

编辑：添加了 MATLAB 代码，以防链接失效。

function FM = fmeasure(Image, Measure, ROI)
%This function measures the relative degree of focus of 
%an image. It may be invoked as:
%
%   FM = fmeasure(Image, Method, ROI)
%
%Where 
%   Image,  is a grayscale image and FM is the computed
%           focus value.
%   Method, is the focus measure algorithm as a string.
%           see 'operators.txt' for a list of focus 
%           measure methods. 
%   ROI,    Image ROI as a rectangle [xo yo width heigth].
%           if an empty argument is passed, the whole
%           image is processed.
%
%  Said Pertuz
%  Abr/2010


if ~isempty(ROI)
    Image = imcrop(Image, ROI);
end

WSize = 15; % Size of local window (only some operators)

switch upper(Measure)
    case 'ACMO' % Absolute Central Moment (Shirvaikar2004)
        if ~isinteger(Image), Image = im2uint8(Image);
        end
        FM = AcMomentum(Image);
                
    case 'BREN' % Brenner's (Santos97)
        [M N] = size(Image);
        DH = Image;
        DV = Image;
        DH(1:M-2,:) = diff(Image,2,1);
        DV(:,1:N-2) = diff(Image,2,2);
        FM = max(DH, DV);        
        FM = FM.^2;
        FM = mean2(FM);
        
    case 'CONT' % Image contrast (Nanda2001)
        ImContrast = inline('sum(abs(x(:)-x(5)))');
        FM = nlfilter(Image, [3 3], ImContrast);
        FM = mean2(FM);
                        
    case 'CURV' % Image Curvature (Helmli2001)
        if ~isinteger(Image), Image = im2uint8(Image);
        end
        M1 = [-1 0 1;-1 0 1;-1 0 1];
        M2 = [1 0 1;1 0 1;1 0 1];
        P0 = imfilter(Image, M1, 'replicate', 'conv')/6;
        P1 = imfilter(Image, M1', 'replicate', 'conv')/6;
        P2 = 3*imfilter(Image, M2, 'replicate', 'conv')/10 ...
            -imfilter(Image, M2', 'replicate', 'conv')/5;
        P3 = -imfilter(Image, M2, 'replicate', 'conv')/5 ...
            +3*imfilter(Image, M2, 'replicate', 'conv')/10;
        FM = abs(P0) + abs(P1) + abs(P2) + abs(P3);
        FM = mean2(FM);
        
    case 'DCTE' % DCT energy ratio (Shen2006)
        FM = nlfilter(Image, [8 8], @DctRatio);
        FM = mean2(FM);
        
    case 'DCTR' % DCT reduced energy ratio (Lee2009)
        FM = nlfilter(Image, [8 8], @ReRatio);
        FM = mean2(FM);
        
    case 'GDER' % Gaussian derivative (Geusebroek2000)        
        N = floor(WSize/2);
        sig = N/2.5;
        [x,y] = meshgrid(-N:N, -N:N);
        G = exp(-(x.^2+y.^2)/(2*sig^2))/(2*pi*sig);
        Gx = -x.*G/(sig^2);Gx = Gx/sum(Gx(:));
        Gy = -y.*G/(sig^2);Gy = Gy/sum(Gy(:));
        Rx = imfilter(double(Image), Gx, 'conv', 'replicate');
        Ry = imfilter(double(Image), Gy, 'conv', 'replicate');
        FM = Rx.^2+Ry.^2;
        FM = mean2(FM);
        
    case 'GLVA' % Graylevel variance (Krotkov86)
        FM = std2(Image);
        
    case 'GLLV' %Graylevel local variance (Pech2000)        
        LVar = stdfilt(Image, ones(WSize,WSize)).^2;
        FM = std2(LVar)^2;
        
    case 'GLVN' % Normalized GLV (Santos97)
        FM = std2(Image)^2/mean2(Image);
        
    case 'GRAE' % Energy of gradient (Subbarao92a)
        Ix = Image;
        Iy = Image;
        Iy(1:end-1,:) = diff(Image, 1, 1);
        Ix(:,1:end-1) = diff(Image, 1, 2);
        FM = Ix.^2 + Iy.^2;
        FM = mean2(FM);
        
    case 'GRAT' % Thresholded gradient (Snatos97)
        Th = 0; %Threshold
        Ix = Image;
        Iy = Image;
        Iy(1:end-1,:) = diff(Image, 1, 1);
        Ix(:,1:end-1) = diff(Image, 1, 2);
        FM = max(abs(Ix), abs(Iy));
        FM(FM<Th)=0;
        FM = sum(FM(:))/sum(sum(FM~=0));
        
    case 'GRAS' % Squared gradient (Eskicioglu95)
        Ix = diff(Image, 1, 2);
        FM = Ix.^2;
        FM = mean2(FM);
        
    case 'HELM' %Helmli's mean method (Helmli2001)        
        MEANF = fspecial('average',[WSize WSize]);
        U = imfilter(Image, MEANF, 'replicate');
        R1 = U./Image;
        R1(Image==0)=1;
        index = (U>Image);
        FM = 1./R1;
        FM(index) = R1(index);
        FM = mean2(FM);
        
    case 'HISE' % Histogram entropy (Krotkov86)
        FM = entropy(Image);
        
    case 'HISR' % Histogram range (Firestone91)
        FM = max(Image(:))-min(Image(:));
        
           
    case 'LAPE' % Energy of laplacian (Subbarao92a)
        LAP = fspecial('laplacian');
        FM = imfilter(Image, LAP, 'replicate', 'conv');
        FM = mean2(FM.^2);
                
    case 'LAPM' % Modified Laplacian (Nayar89)
        M = [-1 2 -1];        
        Lx = imfilter(Image, M, 'replicate', 'conv');
        Ly = imfilter(Image, M', 'replicate', 'conv');
        FM = abs(Lx) + abs(Ly);
        FM = mean2(FM);
        
    case 'LAPV' % Variance of laplacian (Pech2000)
        LAP = fspecial('laplacian');
        ILAP = imfilter(Image, LAP, 'replicate', 'conv');
        FM = std2(ILAP)^2;
        
    case 'LAPD' % Diagonal laplacian (Thelen2009)
        M1 = [-1 2 -1];
        M2 = [0 0 -1;0 2 0;-1 0 0]/sqrt(2);
        M3 = [-1 0 0;0 2 0;0 0 -1]/sqrt(2);
        F1 = imfilter(Image, M1, 'replicate', 'conv');
        F2 = imfilter(Image, M2, 'replicate', 'conv');
        F3 = imfilter(Image, M3, 'replicate', 'conv');
        F4 = imfilter(Image, M1', 'replicate', 'conv');
        FM = abs(F1) + abs(F2) + abs(F3) + abs(F4);
        FM = mean2(FM);
        
    case 'SFIL' %Steerable filters (Minhas2009)
        % Angles = [0 45 90 135 180 225 270 315];
        N = floor(WSize/2);
        sig = N/2.5;
        [x,y] = meshgrid(-N:N, -N:N);
        G = exp(-(x.^2+y.^2)/(2*sig^2))/(2*pi*sig);
        Gx = -x.*G/(sig^2);Gx = Gx/sum(Gx(:));
        Gy = -y.*G/(sig^2);Gy = Gy/sum(Gy(:));
        R(:,:,1) = imfilter(double(Image), Gx, 'conv', 'replicate');
        R(:,:,2) = imfilter(double(Image), Gy, 'conv', 'replicate');
        R(:,:,3) = cosd(45)*R(:,:,1)+sind(45)*R(:,:,2);
        R(:,:,4) = cosd(135)*R(:,:,1)+sind(135)*R(:,:,2);
        R(:,:,5) = cosd(180)*R(:,:,1)+sind(180)*R(:,:,2);
        R(:,:,6) = cosd(225)*R(:,:,1)+sind(225)*R(:,:,2);
        R(:,:,7) = cosd(270)*R(:,:,1)+sind(270)*R(:,:,2);
        R(:,:,7) = cosd(315)*R(:,:,1)+sind(315)*R(:,:,2);
        FM = max(R,[],3);
        FM = mean2(FM);
        
    case 'SFRQ' % Spatial frequency (Eskicioglu95)
        Ix = Image;
        Iy = Image;
        Ix(:,1:end-1) = diff(Image, 1, 2);
        Iy(1:end-1,:) = diff(Image, 1, 1);
        FM = mean2(sqrt(double(Iy.^2+Ix.^2)));
        
    case 'TENG'% Tenengrad (Krotkov86)
        Sx = fspecial('sobel');
        Gx = imfilter(double(Image), Sx, 'replicate', 'conv');
        Gy = imfilter(double(Image), Sx', 'replicate', 'conv');
        FM = Gx.^2 + Gy.^2;
        FM = mean2(FM);
        
    case 'TENV' % Tenengrad variance (Pech2000)
        Sx = fspecial('sobel');
        Gx = imfilter(double(Image), Sx, 'replicate', 'conv');
        Gy = imfilter(double(Image), Sx', 'replicate', 'conv');
        G = Gx.^2 + Gy.^2;
        FM = std2(G)^2;
        
    case 'VOLA' % Vollath's correlation (Santos97)
        Image = double(Image);
        I1 = Image; I1(1:end-1,:) = Image(2:end,:);
        I2 = Image; I2(1:end-2,:) = Image(3:end,:);
        Image = Image.*(I1-I2);
        FM = mean2(Image);
        
    case 'WAVS' %Sum of Wavelet coeffs (Yang2003)
        [C,S] = wavedec2(Image, 1, 'db6');
        H = wrcoef2('h', C, S, 'db6', 1);   
        V = wrcoef2('v', C, S, 'db6', 1);   
        D = wrcoef2('d', C, S, 'db6', 1);   
        FM = abs(H) + abs(V) + abs(D);
        FM = mean2(FM);
        
    case 'WAVV' %Variance of  Wav...(Yang2003)
        [C,S] = wavedec2(Image, 1, 'db6');
        H = abs(wrcoef2('h', C, S, 'db6', 1));
        V = abs(wrcoef2('v', C, S, 'db6', 1));
        D = abs(wrcoef2('d', C, S, 'db6', 1));
        FM = std2(H)^2+std2(V)+std2(D);
        
    case 'WAVR'
        [C,S] = wavedec2(Image, 3, 'db6');
        H = abs(wrcoef2('h', C, S, 'db6', 1));   
        V = abs(wrcoef2('v', C, S, 'db6', 1));   
        D = abs(wrcoef2('d', C, S, 'db6', 1)); 
        A1 = abs(wrcoef2('a', C, S, 'db6', 1));
        A2 = abs(wrcoef2('a', C, S, 'db6', 2));
        A3 = abs(wrcoef2('a', C, S, 'db6', 3));
        A = A1 + A2 + A3;
        WH = H.^2 + V.^2 + D.^2;
        WH = mean2(WH);
        WL = mean2(A);
        FM = WH/WL;
    otherwise
        error('Unknown measure %s',upper(Measure))
end
 end
%************************************************************************
function fm = AcMomentum(Image)
[M N] = size(Image);
Hist = imhist(Image)/(M*N);
Hist = abs((0:255)-255*mean2(Image))'.*Hist;
fm = sum(Hist);
end

%******************************************************************
function fm = DctRatio(M)
MT = dct2(M).^2;
fm = (sum(MT(:))-MT(1,1))/MT(1,1);
end

%************************************************************************
function fm = ReRatio(M)
M = dct2(M);
fm = (M(1,2)^2+M(1,3)^2+M(2,1)^2+M(2,2)^2+M(3,1)^2)/(M(1,1)^2);
end
%******************************************************************

OpenCV 版本的一些示例：

// OpenCV port of 'LAPM' algorithm (Nayar89)
double modifiedLaplacian(const cv::Mat& src)
{
    cv::Mat M = (Mat_<double>(3, 1) << -1, 2, -1);
    cv::Mat G = cv::getGaussianKernel(3, -1, CV_64F);

    cv::Mat Lx;
    cv::sepFilter2D(src, Lx, CV_64F, M, G);

    cv::Mat Ly;
    cv::sepFilter2D(src, Ly, CV_64F, G, M);

    cv::Mat FM = cv::abs(Lx) + cv::abs(Ly);

    double focusMeasure = cv::mean(FM).val[0];
    return focusMeasure;
}

// OpenCV port of 'LAPV' algorithm (Pech2000)
double varianceOfLaplacian(const cv::Mat& src)
{
    cv::Mat lap;
    cv::Laplacian(src, lap, CV_64F);

    cv::Scalar mu, sigma;
    cv::meanStdDev(lap, mu, sigma);

    double focusMeasure = sigma.val[0]*sigma.val[0];
    return focusMeasure;
}

// OpenCV port of 'TENG' algorithm (Krotkov86)
double tenengrad(const cv::Mat& src, int ksize)
{
    cv::Mat Gx, Gy;
    cv::Sobel(src, Gx, CV_64F, 1, 0, ksize);
    cv::Sobel(src, Gy, CV_64F, 0, 1, ksize);

    cv::Mat FM = Gx.mul(Gx) + Gy.mul(Gy);

    double focusMeasure = cv::mean(FM).val[0];
    return focusMeasure;
}

// OpenCV port of 'GLVN' algorithm (Santos97)
double normalizedGraylevelVariance(const cv::Mat& src)
{
    cv::Scalar mu, sigma;
    cv::meanStdDev(src, mu, sigma);

    double focusMeasure = (sigma.val[0]*sigma.val[0]) / mu.val[0];
    return focusMeasure;
}

不保证这些措施是否是解决您的问题的最佳选择，但是如果您找到与这些措施相关的论文，它们可能会给您带来更多见解。希望您发现该代码有用！我知道我做到了。

Answer 4

以耐克的答案为基础。使用 opencv 实现基于拉普拉斯的方法很简单：

short GetSharpness(char* data, unsigned int width, unsigned int height)
{
    // assumes that your image is already in planner yuv or 8 bit greyscale
    IplImage* in = cvCreateImage(cvSize(width,height),IPL_DEPTH_8U,1);
    IplImage* out = cvCreateImage(cvSize(width,height),IPL_DEPTH_16S,1);
    memcpy(in->imageData,data,width*height);

    // aperture size of 1 corresponds to the correct matrix
    cvLaplace(in, out, 1);

    short maxLap = -32767;
    short* imgData = (short*)out->imageData;
    for(int i =0;i<(out->imageSize/2);i++)
    {
        if(imgData[i] > maxLap) maxLap = imgData[i];
    }

    cvReleaseImage(&in);
    cvReleaseImage(&out);
    return maxLap;
}

将返回一个短路，指示检测到的最大清晰度，这是基于我对现实世界样本的测试，这是一个很好的指示相机是否对焦的指标。毫不奇怪，正常值与场景相关，但比 FFT 方法要少得多，FFT 方法必须具有很高的误报率才能在我的应用程序中有用。

Answer 5

我想出了一个完全不同的解决方案。
我需要分析视频静止帧以找到每 (X) 帧中最清晰的帧。这样，我就可以检测运动模糊和/或失焦图像。

我最终使用了 Canny Edge 检测，并且对几乎所有类型的视频都获得了非常非常好的结果（使用 nikie 的方法，我在数字化 VHS 视频和重交错视频方面遇到了问题）。

我通过在原始图像上设置感兴趣区域 (ROI) 来优化性能。

使用 EmguCV ：

//Convert image using Canny
using (Image<Gray, byte> imgCanny = imgOrig.Canny(225, 175))
{
    //Count the number of pixel representing an edge
    int nCountCanny = imgCanny.CountNonzero()[0];

    //Compute a sharpness grade:
    //< 1.5 = blurred, in movement
    //de 1.5 à 6 = acceptable
    //> 6 =stable, sharp
    double dSharpness = (nCountCanny * 1000.0 / (imgCanny.Cols * imgCanny.Rows));
}

Answer 6

感谢 nikie 提出的拉普拉斯建议。
OpenCV 文档 向我指出了相同的方向：
使用 python、cv2 (opencv 2.4.10) 和 numpy...

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
numpy.max(cv2.convertScaleAbs(cv2.Laplacian(gray, 3)))

结果在 0-255 之间。我发现任何超过 200 左右的东西都非常清晰，而到了 100 左右，它就明显模糊了。即使完全模糊，最大值也不会真正低于 20。

Answer 7

我目前使用的一种方法是测量图像中边缘的扩散。寻找这篇论文：

@ARTICLE{Marziliano04perceptualblur,
    author = {Pina Marziliano and Frederic Dufaux and Stefan Winkler and Touradj Ebrahimi},
    title = {Perceptual blur and ringing metrics: Application to JPEG2000,” Signal Process},
    journal = {Image Commun},
    year = {2004},
    pages = {163--172} }

它通常是付费墙后面的，但我见过一些免费的副本。基本上，他们定位图像中的垂直边缘，然后测量这些边缘的宽度。平均宽度给出图像的最终模糊估计结果。较宽的边缘对应于模糊的图像，反之亦然。

这个问题属于无参考图像质量估计领域。如果你在谷歌学术上查找，你会得到很多有用的参考资料。

编辑

这是 nikie 帖子中 5 张图像获得的模糊估计图。值越高，模糊程度越高。我使用固定大小的 11x11 高斯滤波器并改变标准差（使用 imagemagick 的 convert 命令来获取模糊图像）。

如果您比较不同尺寸的图像，请不要忘记按图像宽度进行标准化，因为较大的图像将具有更宽的边缘。

最后，一个重要的问题是区分艺术模糊和不需要的模糊（由焦点缺失、压缩、拍摄对象与相机的相对运动引起），但这超出了像这样的简单方法的范围。举个艺术模糊的例子，看看 Lenna 的图像：Lenna 在镜子中的倒影是模糊的，但她的脸却完全清晰。这有助于对 Lenna 图像进行更高的模糊估计。

Answer 8

上面的答案阐明了很多事情，但我认为进行概念上的区分是有用的。

如果您拍摄一张模糊图像的完美对焦图片会怎样？

只有当您有参考时，模糊检测问题才能很好地解决。参考。例如，如果您需要设计自动对焦系统，您可以比较以不同模糊或平滑程度拍摄的一系列图像，并尝试找到该组中模糊最小的点。换句话说，您需要使用上面说明的技术之一交叉引用各种图像（基本上 - 在方法中进行各种可能的细化级别 - 寻找具有最高高频内容的一张图像）。

Answer 9

我尝试了基于 this 帖子中的拉普拉斯滤波器的解决方案。这对我没有帮助。因此，我尝试了 这篇 帖子中的解决方案，这对我的情况很有好处（但很慢）：

import cv2

image = cv2.imread("test.jpeg")
height, width = image.shape[:2]
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

def px(x, y):
    return int(gray[y, x])

sum = 0
for x in range(width-1):
    for y in range(height):
        sum += abs(px(x, y) - px(x+1, y))

模糊程度较低的图像具有最大 sum 值！

您还可以通过更改步骤来调整速度和准确性，例如，

这部分

for x in range(width - 1):

您可以用这个替换

for x in range(0, width - 1, 10):

Answer 10

已在备受推崇的期刊（IEEE Transactions on Image Processing）上发表的两种方法的 Matlab 代码可在此处获取：https:// /ivulab.asu.edu/software

检查 CPBDM 和 JNBM 算法。如果你检查代码，移植并不难，顺便说一下，它是基于 Marzialiano 的方法作为基本功能的。

Answer 11

我在matlab中使用fft实现了它，并检查fft计算平均值和std的直方图，但也可以完成拟合函数

fa =  abs(fftshift(fft(sharp_img)));
fb = abs(fftshift(fft(blured_img)));

f1=20*log10(0.001+fa);
f2=20*log10(0.001+fb);

figure,imagesc(f1);title('org')
figure,imagesc(f2);title('blur')

figure,hist(f1(:),100);title('org')
figure,hist(f2(:),100);title('blur')

mf1=mean(f1(:));
mf2=mean(f2(:));

mfd1=median(f1(:));
mfd2=median(f2(:));

sf1=std(f1(:));
sf2=std(f2(:));

Answer 12

这就是我在 Opencv 中检测某个区域的焦点质量所做的事情：

Mat grad;
int scale = 1;
int delta = 0;
int ddepth = CV_8U;
Mat grad_x, grad_y;
Mat abs_grad_x, abs_grad_y;
/// Gradient X
Sobel(matFromSensor, grad_x, ddepth, 1, 0, 3, scale, delta, BORDER_DEFAULT);
/// Gradient Y
Sobel(matFromSensor, grad_y, ddepth, 0, 1, 3, scale, delta, BORDER_DEFAULT);
convertScaleAbs(grad_x, abs_grad_x);
convertScaleAbs(grad_y, abs_grad_y);
addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0, grad);
cv::Scalar mu, sigma;
cv::meanStdDev(grad, /* mean */ mu, /*stdev*/ sigma);
focusMeasure = mu.val[0] * mu.val[0];