Nik*_*iki 151
估算图像清晰度的另一种非常简单的方法是使用拉普拉斯(或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对3.3),主要是因为这种方法仅使用图像中最强的对比度,而FFT基本上是整个图像的平均值.功能也在下降得更快,因此设置"模糊"阈值可能更容易.
- +1.你能为你的情节提供一个传奇吗？:P (52认同)
- @Drazick:我不确定这是否可能.例如,看看Lena图像:虽然区域是焦点,但是有很大区域没有对比(例如Lena的皮肤).我想不出一种方法来判断这样一个平滑的区域是否"模糊",或者将它与一个离焦区域区分开来.您应该将此问题作为一个单独的问题(可能在DSP.SE上).也许其他人有更好的想法. (4认同)
- 它适合运动模糊吗？还是只用于像高斯那样的模糊？ (2认同)
Sim*_*got 126
是的.计算快速傅立叶变换并分析结果.傅里叶变换可以告诉您图像中存在哪些频率.如果存在少量高频,则图像模糊.
定义术语"低"和"高"取决于您.
编辑:
如评论中所述,如果您想要一个表示给定图像模糊度的浮点数,则必须计算出合适的指标.
nikie的答案提供了这样一个指标.使用拉普拉斯内核卷积图像:
1
1 -4 1
1
并在输出上使用强大的最大度量标准来获取可用于阈值处理的数字.尝试在计算拉普拉斯算子之前避免过多地平滑图像,因为你只会发现平滑的图像确实模糊:-).
- 唯一的问题是'低'和'高'也是场景相关的.+1 (8认同)
- 只有一点是非常重要的,你必须知道(至少粗略地)你预期的*预模糊*图像(频率)内容是什么.这是正确的,因为频谱将是原始图像的频谱乘以模糊滤波器的频谱.因此,如果原始图像已经主要具有低频率,那么如何判断它是否模糊？ (6认同)
- 除非您的图像是循环的,否则图像边界通常会有锐边,导致频率非常高 (4认同)
- 您通常会虚拟地扩展图像以避免这种影响.你也可以用小窗户计算当地的fft. (2认同)
- 如果您拍摄一张空白的白色图表，您将无法判断图像是否模糊。我认为 OP 需要一些绝对的锐度测量。预模糊图像可能根本不存在。您必须花点功夫才能获得正确的指标，但 fft 可以帮助解决此问题。从这个角度来看，nickie 的回答比我的好。 (2认同)
mev*_*ron 76
在使用自动对焦镜头的一些工作中,我遇到了这组用于检测图像焦点的非常有用的算法.它在MATLAB中实现,但大多数功能很容易通过filter2D移植到OpenCV .
它基本上是许多焦点测量算法的调查实现.如果您想阅读原始论文,请在代码中提供对算法作者的参考.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;
}
无法保证这些措施是否是您问题的最佳选择,但如果您追查与这些措施相关的文件,它们可能会为您提供更多见解.希望您找到有用的代码!我知道我做到了.
Yau*_*aur 31
建立耐克的答案.它直接用opencv实现基于laplacian的方法:
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方法要小得多,因为在我的应用中,误差率很高.
Gol*_*k84 21
我提出了一个完全不同的解决方案.我需要分析视频静止帧以找到每个(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));
}
- 这效果惊人（我在这个问题中测试过的最好的一个）。这是 python 实现，供有需要的人使用。`image_canny = cv2.Canny(image, 175, 225);nonzero_ratio = numpy.count_nonzero(image_canny) * 1000.0 / image_canny.size;返回nonzero_ratio;` (3认同)
gge*_*z44 16
感谢nikie对拉普拉斯的伟大建议. OpenCV文档向我指出了同样的方向:使用python,cv2(opencv 2.4.10)和numpy ......
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
numpy.max(cv2.convertScaleAbs(cv2.Laplacian(gray_image,3)))
结果在0-255之间.我发现任何超过200的东西都非常关注,到100时,它显然是模糊的.即使它完全模糊,最大值也永远不会低于20.
- 我的3张照片中有255张.对于一张完美聚焦的照片,我得到108.所以,我认为该方法的有效性取决于某些东西. (3认同)
我目前使用的一种方法是测量图像中边缘的扩散.寻找这篇论文:
@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} }
它通常在付费墙后面,但我看到了一些免费副本.基本上,他们在图像中定位垂直边缘,然后测量这些边缘的宽度.平均宽度给出图像的最终模糊估计结果.较宽的边缘对应于模糊的图像,反之亦然.
该问题属于无参考图像质量估计的领域.如果您在Google学术搜索中查找,您将获得大量有用的参考资料.
编辑
这是nikie帖子中5幅图像的模糊估计图.较高的值对应于较大的模糊.我使用固定大小的11x11高斯滤波器并改变标准偏差(使用imagemagick的convert命令来获得模糊图像).
如果您比较不同尺寸的图像,请不要忘记按图像宽度进行标准化,因为较大的图像将具有较宽的边缘.
最后,一个重要的问题是区分艺术模糊和不期望的模糊(由焦点未命中,压缩,主体与相机的相对运动引起),但这超出了像这样的简单方法.作为一个艺术模糊的例子,看看Lenna的形象:Lenna在镜子里的反射是模糊的,但她的脸完全是焦点.这有助于Lenna图像的更高模糊估计.
上面的答案阐明了很多事情，但我认为进行概念上的区分是有用的。
如果您拍摄一张模糊图像的完美对焦照片会怎样？
只有当您有参考时，模糊检测问题才能很好地提出。例如，如果您需要设计自动对焦系统，您可以比较以不同模糊或平滑程度拍摄的一系列图像，并尝试找到该组中模糊最小的点。换句话说，您需要使用上面说明的技术之一交叉引用各种图像（基本上 - 在方法中进行各种可能的细化级别 - 寻找具有最高高频内容的一张图像）。
- 换句话说，这是一个相对的概念，只能判断一张图像是否比另一张类似图像更加模糊或更少模糊。即，其 FFT 中是否具有或多或少的高频内容。特殊情况：如果图像的相邻像素具有最大和最小亮度怎么办？例如，全黑像素紧邻全白像素。在这种情况下，它是一个完美的焦点，否则从黑色到白色的过渡会更平滑。摄影中不太可能实现完美对焦，但问题没有指定图像的来源（可能是计算机生成的）。 (2认同)
我从这篇文章中尝试了基于拉普拉斯滤波器的解决方案。它没有帮助我。所以，我尝试了这篇文章中的解决方案，这对我的情况有好处（但速度很慢）：
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):