Efficient blur estimation using multi-scale quadrature filters

Signal Processing 93 (2013) 1988–2002

Contents lists available at SciVerse ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Efficient blur estimation using multi-scale quadrature filters Seyfollah Soleimani a,b,n, Filip Rooms a, Wilfried Philips a,1 a b

TELIN-IPI-iMinds, Ghent University, St-Pietersnieuwstraat 41, B-9000 Gent, Belgium Department of Computer Engineering, Faculty of Engineering, Arak University, Arak 38156-8-8349, Iran

a r t i c l e in f o

abstract

Article history: Received 12 April 2012 Received in revised form 27 October 2012 Accepted 21 December 2012 Available online 1 January 2013

Blur estimation is required in image processing techniques such as auto-focussing, quality assessment for compressed images and image fusion. In this paper a multi-scale local blur estimation method is proposed. We use the energy of a pair of quadrature filters with first and second derivatives of a Gaussian at several scales as its constituents. A new strategy for analyzing the extrema of energy across scale is proposed. Comparing to the methods which use just a Gaussian first derivative kernel, a smaller number of scales needs to be processed. Also our method yields only one response at the centroid of line-shape features at a position independent of the scale. This is in contrast to other methods which yield multiple responses at scale dependent positions. We evaluated the method for synthetic and real images from the LIVE database. Depending on details of the image, the proposed method is several to tens of times faster in comparison with using just a Gaussian first derivative. The accuracy of the blur estimation is found to be the best or second best while comparing with 16 different methods for Gaussian blur. In addition, the performance is still good for motion blurred images. & 2012 Elsevier B.V. All rights reserved.

Keywords: Multi-scale edge detection Blur estimation Quadrature filters Wavelet transform

1. Introduction Blur estimation is required in some image processing techniques such as auto-focusing, quality assessment of compression methods and image fusion [1–3]. This technique has been explored by various authors in the past. It is possible to categorize these methods in several ways. One of the ways is based on whether or not they use a reference image. In this respect methods are divided into full-reference, no-reference and reduced reference [4]. In full-reference methods the processed image is compared with a reference image. Whereas in no-reference methods, n Corresponding author at: TELIN-IPI-iMinds, Ghent University, St-Pietersnieuwstraat 41, B-9000 Gent, Belgium. Tel.: þ32 9 264 34 12; fax: þ32 9 264 42 95. E-mail addresses: [email protected], [email protected] (S. Soleimani), [email protected] (F. Rooms), [email protected] (W. Philips). 1 Senior member of IEEE.

0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2012.12.015

only the single input image is processed. With reduced reference methods, there is no-reference image, however, some features of the reference image are assumed to be known. For example, when an image is transmitted over a network, it is possible to extract some simple features from the image and send them along with the compressed image. Another way of categorizing blur estimation methods is according to whether they compute a single estimate for whole image or local estimates. The former is known as global methods, while the latter is known as local. Yet, another way of categorizing is methods that operate on blurred images of different scenes versus methods that operate on multiple images of the same scene. It is important that blur estimation algorithms are robust to noise. Methods that use edge regions of images are expected to be more robust against noise than the methods using global characteristics indiscriminately. Some authors have proposed to use multi-scale processing for edge detection and blur estimation. Multi-scale processing enables us to work in several scales instead of

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

just one scale. Coarse scales are less sensitive to noise while they eliminate fine structures and vice versa for fine scales. Following Mallat [5], who proposed to estimate Lipschitz exponents of singularities using the evolution of local maxima of wavelet coefficients across scales, Ducottet proposed a method for edge detection and blur estimation [6]. His method is based on modelling three types of singularities in images: transitions, lines and peaks. These three singularities are modelled as a Heaviside step edge, a one-dimensional delta point feature and a two-dimensional delta line feature, each smoothed by a Gaussian. In his method, first the wavelet transform of the models is computed using the first derivative of the Gaussian as the wavelet. Then by analyzing the local maxima of coefficients in wavelet domain, for each model a closed form formula is derived which is called the maxima function. In other words, every maxima function is created by connecting local maxima of wavelet coefficients associated with a particular point ðx0 ,y0 Þ across scale. The scale is the standard deviation of the Gaussian wavelet (which is denoted by s and determines the width of the wavelet) and ranges from smin to smax with steps of Ds (Table 1 lists notations and definitions used in this paper). Having an analytical maxima function for every model, Ducottet’s method performs as follows: 1. The wavelet transform of the input image is calculated at several scales. 2. The local maxima of wavelet coefficients are located in every scale. 3. Starting from the finest scale, the corresponding local maxima of wavelet coefficients from different scales are combined into maxima functions. Table 1 Denotations and definitions. Symbol

Definition

s

Scale ¼the standard deviation of the Gaussian which is used for creating the scale space Gs ðx,yÞ A Gaussian function with standard deviation of s which is used in creating the scale space Ns Number of processing scales Smin Minimum (finest) processing scale Smax Maximum (coarsest) processing scale Ds The scale step s Blur level of the singularities Gs ðx,yÞ A Gaussian function with standard deviation of s which is used in modelling of the singularities HðxÞ ¼ Hðx,yÞ The Heaviside function 1D Dirac function dðxÞ T s ðx,yÞ Transition model with blur level of s Ls ðx,yÞ Line model with blur level of s Rx 2 erf(x) p2ffiffiffi 0 expðt Þ dt p

n

f ðx,yÞ Nf

smax ET s ðx,y,sÞ ELs ðx,y,sÞ MET s ðsÞ MELs ðsÞ STD MAPE

Convolution operator A two-dimensional function (image) The minimum number of values that an extrema function should contain to keep for processing Maximum blur level in the image Energy of the transition model at scale s Energy of the line model at scale s Extrema function of the transition model Extrema function of the line model Standard deviation of noise or blur Mean absolute percentage error

1989

4. The best fit of analytical maxima functions to every extracted maxima function from input image specifies the type and the blur level associated with the extracted maxima function. 5. The average blurriness of all detected edges is calculated to have a number for whole image.

Ducottet’s method performs well but it is rather slow which is due to use of the Gaussian first derivative which causes two problems: (1) for some features in the image, it produces more than one local maximum in the wavelet domain and (2) the locations of those local maxima shift with scale. So, for a given local maximum of the wavelet coefficients in the current scale, finding its corresponding local maximum in the next coarser scale requires a time consuming neighborhood search. To be able to limit the search to only the eight neighboring pixels, Ducottet imposed Ds r 0:5. Thus, a large number of scales need to be processed. In addition, because the features produce multiple responses, more processing is needed comparing to the situation in which just one response is produced. The occurrence of multiple responses and the high number of scales to be processed cause a high computational cost for the Ducottet method. We used the Ducottet method (with some changes) for blur estimation in our previous work [1], where we reduced the time complexity by setting a scale-dependent threshold on wavelet coefficients in every scale. With the proposed method of this paper, first, the limitation on the scale step is relaxed which reduces the number of processing scales. Second, for any feature, only one response is processed in any scale of the wavelet domain. Therefore, the computational cost diminishes considerably. The experimental results of Section 6 will show that, depending on the details of the image, our method is several to 10-fold faster. In this paper, we explore the multi-scale quadrature filters that use the first and second derivatives of Gaussian filters. Using this type of filter for edge detection along with a special strategy for combining local extrema of filtered signal across the scale, any type of singularities yields one response at the center in any scale. We use the word extrema instead of maxima because in some scales we have to take into account both local minima and local maxima (see Section 3.4). Also, we compare the performance of our method to that of 16 other no-reference blur metrics. The comparison focuses on how much the blur estimates varying monotonically with the true blurring level [10,18]. Any deviation from this rule is considered as an error. We calculate it in such a way that larger deviation is more penalized. For example, suppose that the real blur levels are B1 , B2 , B3 and B1 oB2 o B3 . Also assume that the ^ but estimated blur levels are B^1 , B^2 , B^3 . If B^2 5B3 ^ and ^ ^ ^ B1 o B3 , it is counted as one error and if B2 5B3 B^1 5B^3 , it is counted as two errors. Moreover, we use Spearman rank-order correlation for testing the monotonicity whose results are consistent with the number of errors.

1990

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

The major contributions of this paper can be summarized as follows: 1. Use quadrature filters instead of only first order derivatives of Gaussian for multi-scale edge detection and blur estimation. The employed quadrature filters use both first and second order derivatives of Gaussian. 2. Adjust the algorithm to deal with the new filters. 3. It is a huge improvement of the method proposed in [6] in the sense of speed while the accuracy is better. 4. We compare the performance of the proposed method with that of 16 other no-reference blur metrics for different amounts of noise. 5. For calculating the convolution, Ducottet [6] used a Fourier domain method. However, we used the spatial method of [24], which is faster (the main reason of efficiency is because of using the proposed multi-scale quadrature filters). The paper is organized as follows: in Section 2, we review 16 existing no-reference methods that we will use to compare with our method. In Section 3, the main idea of our proposed method is introduced. In Section 4, the algorithm and implementation details are explained. The reasons why our method is faster are summarized in Section 5. The evaluation of the method and discussion are done in Section 6 and we conclude in Section 7.

 Laplacian(LAP) [10]: The Laplacian operator is based on the second derivatives of the image in the x- and y-directions: First the image is filtered by the following mask: 2 3 0 1 0 6 7 4 1 4 1 5 0 1 0 Then, each pixel in the resulting image is squared and finally the average is calculated as the measure. The wellknown drawback of this metric is its high sensitivity to the noise.

2.2. Statistical algorithms

 Variance (VAR) [9]: It is simply the variance of the whole image: Bvar ¼



M X N 1 X ½f ðx,yÞf 2 MN x ¼ 1 y ¼ 1

where f is the mean of intensities in the image. Although this method is the simplest and most cited measure, it has limited discrimination power [18]. Auto-Correlation (AUTO) [10]: The auto-correlation function is Cðm,nÞ ¼

Mm X X Nn

f ðx,yÞf ðx þm,y þ nÞ

x¼1y¼1

2. Review of other no-reference methods To estimate the amount of blurring of an image, if there is a reference (which is usually the ground truth image), the blurred image can be compared either by extracted features from both images or by directly comparing the gray values. But when there is no reference image, the decision should be made just based on the input image(s). Most methods estimate the blurriness by quantifying the low contrast (low frequencies) or high contrast (high frequencies) regions in the image. A good method should produce an estimate that varies monotonically with respect to the blurring level and also should be robust against noise. In this section, 16 existing no-reference blur metrics are reviewed. We will compare our method with these methods in Section 6. Here, we also introduce an acronym for every method (shown in parentheses) to simplify referring to them. We categorized them roughly into four groups:







2.1. Derivative-based algorithms

 The Normalized Gradient (GRA) [10]: This metric is based on the first order derivative. If the image is f ðx,yÞ with size of M  N, then the metric is Bgrad,n ¼

X qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X Nn 1 Mn ½f ðx þ n,yÞf ðx,yÞ2 þ ½f ðx,yþ nÞf ðx,yÞ2 MN x ¼ 1 y ¼ 1

In [10], n proposed to be a number less than 10. This metric inherits the sensitivity of the derivation operation to the noise.



where 1 r mr M and 1 rn rN. In this method, some values for m and n are selected and then computed auto-correlation coefficient is divided by the sum of squares of all intensities in the image. In our experiments m and n were set to 3. Natural Scene Statistics for JPEG2K (NSS) [7]: This method uses natural scene statistics for quality assessment based on the reasoning that those statistics and image quality are related. So, they estimate the metric as the change in natural behavior of the image. Then the metric is calibrated against subjective scores. The kurtosis based metric by Zhang (ZHANG) [12]: Kurtosis is defined as a descriptor of the shape of a probability distribution function. A wider distribution has smaller kurtosis. In this method, the spectral density function of the image is considered as a probability density function and its kurtosis is used as a sharpness metric. Blurry images produce larger kurtosis. The Cumulative Probability Blur Detection (CPBD) [17]: In this method, first edge detection is performed and then for every detected edge, the probability of blur is estimated. Afterwards, a probability density function (Pdf) is derived from Pdfs of detected edges. Finally the cumulative probability of blur is calculated from the derived Pdf. The Just Noticeable Blur (JNB) metric [4]: A perceptual sharpness metric is derived based on measured justnoticeable blurs and probability summation over space, which takes into account the response from HVS (human vision system) to sharpness at different contrast levels. This method gives the relative amount of blurriness in images with different contents.

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

2.3. Transform-domain algorithms

1991

2.4. Other algorithms

 Spectrum [2]: This metric has been used for auto-

 No-reference JPEG (NRJPEG) [8]: In this method,

focussing. First the Fourier transform of the image is calculated. Then the DC component of the image is set to zero (if the input image is f ðx,yÞ by size M  N and its discrete Fourier transform is Fðu,vÞ by the same size, by considering the lowest frequencies at ðu,vÞ ¼ ð1,1Þ, the DC component of the image is Fð1,1Þ). The DC component of the image is the average of all intensities in the image. It is assumed that by changing the focus, the DC component of the image does not change significantly, so it is ignored. When an image comes into focus, low frequency components decrease and when the image becomes blurry, the low frequency components increase. In this method the summation of a part of low frequency components is calculated as the metric. In our experiments, the metric is estimated by the summation of frequency components in the first quarter. It means that by increasing the sharpness, this metric reduces. This metric is sensitive to noise as well because it takes into account all information in the image. Marichal [11]: This method is based on the histogram of non-zero DCT coefficients in 8  8 blocks of the image. The idea originates from the fact that DCT coefficients reflect the frequency distribution of an image block and when an image is more blurred, produces more zero coefficients. Ratio of High-pass and Low-pass Bands of Wavelet Transform (RHLWT) [18]: If the input image is f ðx,yÞ with size of M  N, in this method, first, the discrete wavelet transform of the image is calculated using an orthogonal wavelet (like Daubechies wavelet). This transform produces a low-pass band (lw(f)) and several high-pass bands (hw(f)). Then the sharpness metric is calculated as follows:

two measures are combined to estimate the blurriness of the image. The measures are calculated horizontally and vertically. The first measure is the blockiness which is estimated as the average differences across block boundaries. The horizontal blockiness is





1 Jhw ðf ÞJ B ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 d m 1 Jlw ðf ÞJ







where d is the depth of decomposition and m is the dilation factor of the wavelet and J  J denotes the discrete Euclidean norm. The Riemannian tensor based metric (Riemannian) [16]: It is based on the idea of viewing gray level images as surfaces in a non-Euclidean 3D space (and color images in a 5D space) which comes from a Riemannian differential geometry. The derived sharpness metric is the average of determinant tensor metric of all pixels in the image. The tensor metric is a 2  2 matrix whose elements are some combinations of image first derivatives in the x- and y-directions. The Kurtosis of Wavelet Coefficients (Wavelet) [15]: This method is based on measuring the kurtosis in the wavelet domain. First the dyadic wavelet transform of the image is calculated and then the kurtosis is applied on the vertical and horizontal detail sub-bands. Then, the average of two computed metrics is calculated as the final metric. High-pass-to-band-pass frequency ratio (HP) [14]: The metric is based on the high-pass-to-band-pass frequency ratio applied to local features that are extracted by thresholding the bandpass filter output.

Bh ¼

M bN=8c1 X X 1 9dh ði,8jÞ9 MðbN=8cÞ i ¼ 1 j ¼ 1

where dh ðm,nÞ ¼ f ðm,n þ 1Þf ðm,nÞ. The second measure is called the activity measure which in turn is calculated using two other factors. The first factor is the average absolute difference between in-block image samples 2 3 M N1 X X 14 8 Ah ¼ 9d ði,jÞBh 95 7 MðN1Þ i ¼ 1 j ¼ 1 h and the second factor is the zero-crossing rate which is estimated by Zh ¼

M N 2 X X 1 z ðm,nÞ MðN2Þ i ¼ 1 j ¼ 1 h

while zh ðm,nÞ is 1 for horizontal zero-crossing at dh ðm,nÞ and is 0 otherwise. The same features are calculated vertically which are denoted by Bv, Av and Zv. Then the overall features are computed by B¼

Bh þ Bv , 2

A¼

Ah þ Av , 2

Z¼

Zh þ Zv 2

Finally the above features are combined to reach the proposed metric: S ¼ a þ bBg1 Ag2 Z g3



where a, b, g1 , g2 and g3 are estimated with the subjective tests. Marziliano Metric [13]: In this method, first edge detection is performed and then the width of any detected edge is estimated. The width of any detected edge is considered as the distance between the two extrema on both sides of the edge pixel. Then the average of widths of all edge pixels is calculated as the sharpness metric.

3. Energy of quadrature filters In this section we elaborate the proposed method for edge detection and blur estimation. First let us model singularities of images.

3.1. Singularities models We model singularities of images just as transitions (T s ) and lines (Ls ) which are defined as follows:    A x pffiffiffi 1 þerf T s ðx,yÞ ¼ AHðxÞnGs ðx,yÞ ¼ 2 s 2 Ls ðx,yÞ ¼ 2ps2 AdðxÞnGs ðx,yÞ ¼ 2ps2 AGs ðx,0Þ

1992

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

where n is the convolution, A is the amplitude, H is the Heaviside function: ( 0 if x r0 HðxÞ ¼ 1 if x 40 and Gs ðx,yÞ ¼

 2  x þ y2 exp  2 2 2s 2ps

neighborhood in any scale increases and if the scale step reduces, the number of scales should be increased to span the same scale range. In other words, using linear filters, it is not possible to reduce the computational cost.

3.3. Evolution of energy of quadrature filters over scales

1

The erf function is defined by the following integral: Z x 2 erfðxÞ ¼ pffiffiffiffi expðt 2 Þ dt

p

The energy of quadrature filters [19–22] is a possible solution to overcome the drawback of shifting problem in

0

The models are normalized to have peak amplitude A at (0,0). In our discussion s (standard deviation of the Gaussian) is the unknown blur level of the singularities (it should not be confused with s, which is the standard deviation of a Gaussian convolving with the whole image to produce scale space).

3.2. Linear filters and shifting problem across scale For edge detection, if we convolve the input image with the Gaussian first derivative filter (which is an odd filter) in the x- and y-directions and calculate the modulus of the responses, the line-type features (which are even features) will produce two parallel lines of local maxima located on each side of the line. The locations of these local maxima depend on the scale (s) of the filter. This problem is illustrated in Fig. 1 in one dimension and in Fig. 2 in two dimensions. In Fig. 1(a) the profile of a line-type feature is shown and in Fig. 1(b) the local maxima of modulus of filtered signal in two different scales are shown (one pair by dash curve and the other by solid curve). In Fig. 2(a), a line-type singularity is depicted and its local maxima of filtered signal in two scales are shown in Fig. 2(b) (the local maxima of one scale by ‘þ’ symbol and in the other scale by ‘O’ symbol). From these two figures, it is clear that the responses for every scale are at different locations. If we use the Gaussian second derivative filter, this problem does not occur but a similar problem occurs for transition-type features (which are odd features). So, using linear filters for edge detection, to find corresponding local maxima across scale, we have to search a neighborhood whose size is proportional to the scale step. Hence, if the scale step increases, the searching

Fig. 2. (a) A line-type edge (any point presents a gray value). (b) Local modulus maxima of wavelet coefficients of (a) using the Gaussian first derivative as the wavelet in two different scales. One group is shown by circles and the other group by plus symbols.

Fig. 1. (a) The profile of a line-type edge. (b) Local modulus maxima of wavelet coefficients of (a) using the Gaussian first derivative as the wavelet in two different scales.

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

linear filters. The main property of using the energy of quadrature filters for edge detection is that one of its responses (local extrema) to both transition-type and line-type features always occurs at the centroid of every feature. While the central response is on extremum, it is not necessarily a local maximum but can also be a local minimum. For line-type features, there are still two other responses on both sides of the central response which should be eliminated. Kube and Perona stated that there are two groups of quadrature filters: derivative pairs and Hilbert pairs [20]. They proved in one dimension that under some reasonable assumptions, Hilbert pairs do not have the causality property but derivative pairs have. The causality property means that no new features (like local extrema or zerocrossings) are created when we convolve the input signal with a filter of coarser scale in comparison with finer scales. In their experiment, a pair of the first and second derivatives of the Gaussian is used and edge points in one-dimensional signal are found as local maxima. Here, we use steerable pair of Gaussian first and second derivatives in several scales for edge detection and blur estimation. The proposed method of this paper enables us to process just a central scale-independent response for any feature (in any scale). This reduces the processing time for two reasons in comparison with using linear filters. First, because the number of processing scales is reduced and second, because, in any scale, the number of extrema is reduced. The number of scales is reduced as the scale step is large while for finding the corresponding extrema of coefficients in a scale, we do not need to search a significant region in other scales. In addition, the number of extrema in any scale is reduced, because just the central response has to be processed for any pixel of the features while using linear filters, for any pixel of some features, several responses have to be analyzed. The energy of the image response to the pair of first and second derivatives of the Gaussian is Eðx,y,s, yÞ ¼ ðf ðx,y,s, yÞÞ2 þ ðf ðx,y,s, yÞÞ2 0

00

ð1Þ

the dominant direction of (1). The dominant direction is defined as the orientation y in which the Eðx,y,s, y) is maximal over a given position and scale. The dominant direction and the magnitude of energy (Eðx,y,s, yÞ) are saved for every pixel. Then the local extrema of energy should be found in every scale. Corresponding energy extrema in different scales are called extrema functions (instead of maxima functions of Ducottet). In the next sub-section we derive a closed form formula for extrema functions of the above models. The analytical extrema functions are fitted to extracted extrema functions of the input image to estimate the associated blur levels. 3.4. Scale-independent responses For every model, we should find the direction along which the magnitude of energy (Eq. (1)) is maximum for every (x,y,s). For the transition model, the dominant direction for all (x,y,s) is y ¼ 0, so the energy of transition model (ET s ) is " #   A2 s2 s2 x2 x2 1 þ ET s ðx,y,sÞ ¼ exp  2pðs2 þ s2 Þ s2 þ s2 ðs2 þ s2 Þ2 ð2Þ The extrema of ET s are on a line at x¼0 and always at its maxima. So, the locations of these local maxima are independent of s (scale) and s (blurriness of model). The extrema function for x ¼0 of the transition model is MET s ðsÞ ¼

0

f ðx,y,s, yÞ ¼ f x ðx,y,sÞ cosðyÞ þ f y ðx,y,sÞ sinðyÞ 00

f ðx,y,s, yÞ ¼ f xx ðx,y,sÞ cos2 ðyÞf xy ðx,y,sÞ sinð2yÞ þf yy ðx,y,sÞ sin2 ðyÞ and f x ðx,y,sÞ ¼ ðf ðx,yÞÞns

@Gs ðx,yÞ , @x

f xx ðx,y,sÞ ¼ ðf ðx,yÞÞns2

@2 Gs ðx,yÞ , @x2

f yy ðx,y,sÞ ¼ ðf ðx,yÞÞns2

@2 Gs ðx,yÞ , @y2

f y ðx,y,sÞ ¼ ðf ðx,yÞÞns

@Gs ðx,yÞ @y

@2 Gs ðx,yÞ @x@y  2  1 x þ y2 Gs ðx,yÞ ¼ exp  2 2 2ps 2s f xy ðx,y,sÞ ¼ ðf ðx,yÞÞns2

where s is the scale (standard deviation of the Gaussian) and y is the direction of derivation. We are interested in

A2 s2 2pðs2 þ s2 Þ

ð3Þ

The scale-independent extrema of the Line model are at x ¼0 and the scale-dependent ones are at sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðs2 þ s2 Þ x¼ 7 ð3s2 s2 Þ 7 ð3s2 s2 Þ2 þ4s2 ðs2 2s2 Þ 2 2s The scale-independent responses of Line model are pffiffiffi maxima if s 4 s= 2, otherwise they are minima. The extrema function for x¼0 of the line model is MELs ¼

where

1993

A2 s4 s2 ðs2 þ s2 Þ3

ð4Þ

The extrema functions of the models are shown in Fig. 3 for s ¼ 3 and A ¼1. The energy of transition type features has one local extremum at the center which is scale-independent but for line type features, in addition to a response at the centroid of the feature, there are two other scaledependent extrema. These scale-dependent extrema should be eliminated because the above-derived extrema functions are just valid for scale-independent extrema. To keep only the scale invariant responses, we exploit that the desired response does not shift with scales, whereas the scale-dependent responses do shift with scale. Therefore, when at the same location in several successive scales, a local extremum exists, we infer that there is a scale-independent response at that location. We define a parameter Nf which is the minimum number of scales for which, an extremum should exist at the same location to

1994

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

Fig. 3. Extrema functions of models.

correspond with a valid extrema function. In our experiments we set N f ¼ 4 (while the number of processing scales is Ns ¼ 4). This means that we just process the extrema functions that contain at least 4 values. Our experiments show that by setting N f ¼ 4 and Ds ¼ 3, the scale-dependent responses are eliminated, while scaleindependent ones survive. 4. Algorithm and implementation details The algorithm is outlined in Fig. 4. Fig. 4. Schematic of the algorithm.

1. fx, fy, fxx, fxy and fyy are calculated in some scales. Ducottet calculated them using a Fourier domain method. We also tested the spatial domain method of [24]. For calculating in Fourier domain, we used the library FFTW from [23]. In Section 6, we show the results using both Fourier domain and spatial domain for Ducottet’s method and our method. The experiments show that the spatial domain method is faster than the Fourier domain method. 2. For every pixel in any scale, the direction along which the energy (Eq. (1)) is maximum is found. Using the derivative of Eq. (1) to find the dominant direction amounts to solve a quartic (order 4) equation. Because finding the roots is complicated, we rather just calculate the energy in eight directions:

y ¼ 0 or p-E ¼ f 2x þ f 2xx y ¼ p=4 or 5p=4-E ¼ 0:5ðf x þ f y Þ2 þ 0:25ðf xx þf yy 2f xy Þ2

y ¼ p=2 or 3p=2-E ¼ f 2y þf 2yy y ¼ 3p=4 or 7p=4 -E ¼ 0:5ðf x f y Þ2 þ0:25ðf xx þ f yy þ 2f xy Þ2 Then the maximum value and its direction are saved for every pixel. 3. The extrema in every scale are found. As we said before the scale-independent extrema of transition model are maxima for every s and s, but for the Line model they

pffiffiffi are maxima when s 4 s= 2, otherwise they are minima. If we assume that the maximum blur level pffiffiffi that exists in the image is smax , then for s o smax = 2 we should findpboth local minima and maxima and ffiffiffi when s4 smax = 2 we should only find local maxima. 4. Starting from the coarsest scale, for every local extremum we find its corresponding extremum in the next finer scale and continue this procedure till the finest scale. It is possible to start from the finest scale like Ducottet and some other authors, but starting from the coarsest scale is faster. Because when we start from the finest scale, more extrema functions are created at first which are later rejected. For ideal edge models, the position of scaleindependent extrema will not shift by changing the scale. But in practice, the shapes of the features in natural images are not exactly the same as those of the theoretical models. Thus, the scale-independent extrema will also move slightly. For faster processing, we assume that this displacement is at most one pixel. So, for finding the corresponding extremum at another scale, first we look at the same location. If we do not find any, then we look at the eight neighboring pixels. Here, we should insist on one important point. The implicit assumption of our algorithm is that features are isolated enough so that their responses do not overlap. However, it is possible that in coarser scales, distinct features overlap. When the features start to

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

overlap, the mentioned statements about scale-independent extrema will not be valid anymore. To reject these extrema functions, there are two strategies. The first one is implicit in the algorithm: sometimes when in a coarse scale, there is a response associated with a combined feature, there is no corresponding response for it in finer scale(s) because in finer scales, those features are not combined. So, those extrema will be ignored automatically. The second strategy is based on the fact that the extrema function of the transition model (3) is strictly increasing as a function of s. It is also possible to normalize line extrema function (Eq. (4)) to be strictly increasing. If we multiply the extrema function of the line model (Eq. (4)) by s2 and then apply cubic root to it, we have   s2 NMELs ðsÞ ¼ ðs2 MELs Þ1=3 ¼ ðA2 s2 Þ1=3 2 2 s þs Then we obtain a strictly increasing function. So, for every extracted extrema function we pursue two hypotheses. First it is a transition type and second it is a line type. In any case it should be strictly increasing either by itself or after the mentioned normalization. If it is not so, we reject it. 5. Extrema functions of the models are fitted to the extracted extrema functions for classification to transition or line and also to estimate the blur level. The fitting is done by converting the model extrema functions to linear functions with respect to s2 as follows: For transition hypotheses: rffiffiffi   s2 2p b ¼ 2 ðs2 þ s2 Þ ¼ aS þ b ) s^ ¼ MET s ðsÞ a A For line hypotheses: rffiffiffi    1=3 s2 1 b 2 2 ^ ¼ ðs þ s Þ ¼ a S þ b ) s ¼ NMELs ðsÞ a A2 s2 where S ¼ s2 . 6. The average of blurriness of all detected edges is calculated.

5. Time efficiency We denote the number of processing scales in Ducottet’s method by N D s and the number of processing scales in our method by NQs . Let us assume that there is a line-type feature in the input image with length of l pixels. Based on the previous discussions, the number of extrema in Ducottet’s method in each scale is 2l while in our method it is l. To keep the connection between corresponding extrema in different scales, a tree structure is needed. For any scale, in Ducottet’s method, 2l nodes are inserted in the tree and the depth of the tree is ND s , while using our method, l nodes are inserted in the tree and the depth of the tree is NQs . Therefore, the ratio of processing times of Ducottet’s method to our method for line type features is Q D Q ð2l  N D s Þ=ðl  N s Þ ¼ 2N s =N s .

1995

Table 2 Comparing time complexity of Ducottet’s method and our method. Smin

Smax

DsD

ND s

DsQ

N Qs

Ratio

Line 1 1 1 1

10 10 17 21

1 0.5 0.2 0.1

10 19 81 201

3 3 4 4

4 4 5 6

5 9.5 32.4 67

Transition 1 1 1 1

10 10 17 21

1 0.5 0.2 0.1

10 19 81 201

3 3 4 4

4 4 5 6

2.5 4.75 16.2 33.5

If there is a transition type feature with length of t, both methods have t extrema in any scale, but each extrema function in Ducottet’s method contains N D s values while in our method, each extrema function contains N Qs . In this case, the ratio of processing times of Ducottet’s Q method to our method is N D s =N s . Every extrema function is analyzed to classify to the line or transition types and also to calculate the blur level. As the number of line type extrema functions in Ducottet’s method is 2l and each contains ND s values, the total number of values of line type features using Ducottet’s Q method is 2l  ND s while for our method it is l  N s . So D Q again the ratio 2N s =Ns holds. With the same reasoning Q the ratio for transition type features is ND s =N s . In Table 2 we show selected values for each parameter of both algorithms and resulting processing time ratios. It is clear from this table that our method is several to tens of times faster. 6. Evaluation To evaluate the proposed method, we use four sets of test images. The first set contains synthetic images and the second set contains real images from Live database [3] which are distorted by Gaussian blur. The third set comprises images that are distorted by motion blur and the fourth set contains images that are compressed with the JPEG2000 standard. From now on, the acronym STD is used to indicate the standard deviation of noise or the blur level. Filtering the input image with Gaussian derivatives as a part of process in both our method and Ducottet’s method can be performed in either the spatial domain or the frequency domain. Ducottet implemented his algorithm in the frequency domain. Here, we implemented the Ducottet method and our method in both spatial and frequency domains. For calculating in the spatial domain, we used the code from [24] and for calculating in the frequency domain we used the FFTW library [23]. In the following we use SpaDu for referring to the Ducottet method using spatial domain calculations and FftDu for referring to the Ducottet method using frequency domain calculations. Also, we will use SpaQF for referring to our proposed method using spatial domain calculations and FftQF for referring to our proposed method using frequency domain calculations.

1996

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

images but not for larger values. So in our experiments for Ducottet’s method we set Ds ¼ 1.

6.1. Synthetic images Every image in this set contains a straight horizontal transition and a straight horizontal line with the STDs of blur levels 1.0–10 pixels with 0.5 as the step. Furthermore, every image was distorted with white Gaussian noise with STDs of 5, 10 and 15. So, we made a data set of 76 images (considering noise free images). The proposed blur estimation algorithm was applied to the synthetic images. For every input image, locations of edge pixels along with their type (transition or line) and their blur levels were calculated. Having the real blur level for these synthetic images, we can calculate the mean absolute percentage error (MAPE) of estimation: MAPE ¼

For our method, we set Ds ¼ 3 and by this setting all features were detected correctly. The mean absolute percentage errors (MAPE) for Ducottet’s method and our method are shown in Table 4. The Ducottet method was tested with the number of scales (Ns) of 4, 5, 6 and 8 and our method was tested for Ns ¼ 4. In the frequency domain our method has less error in all tested situations except in the case of noise STD of 5 when Ns ¼ 8. In the spatial domain when N s Z5 for noise free images, Ducottet’s method has a bit less error. For other amounts of noise, Ducottet’s method is better when Ns Z 6. Anyway the error of our proposed method in the spatial domain is at most 6.6% which is a good estimation. Based on these experiments we chose to use 6 or more scales for Ducottet’s method for the best in the real images in the next subsection. The average blur levels for all detected edges using our method by calculating the frequency domain and spatial domain are plotted in Fig. 5(a) and (b), respectively. As can be seen the proposed method using both domains works well even in the presence of heavy noise for a large range of blurriness. It confirms that our method is intrinsically robust against noise.

F 9B^i Bi 9 100 X F i ¼ 1 Bi

B^i is the estimated blur level of the ith image, Bi is the real blur level of the ith image and F is the number of images. The parameters were set as shown in Table 3. To reduce the effect of noise, we increased the finest scale (smin) by increasing the level of noise because finer scales are more sensitive to noise. To check the best value for scale step in Ducottet’s method, some values were tested:

 By setting Ds ¼ 3, Smin ¼ 1, Ns ¼ 4 and Nf ¼ 4, all 



6.2. Real images

transition-type features were detected but no linetype feature was detected. By setting Ds ¼ 2, Smin ¼ 1, Ns ¼ 4 and N f ¼ 4, all transition-type features were detected but line-type features in the majority of images (10 out of 19) were not detected. By setting Ds ¼ 1, Smin ¼ 1, Ns ¼ 8 and N f ¼ 8, every type of feature was detected. It means that the Ducottet method works well by a scale step of 1 for synthetic

The second set of test images was drawn from LIVE database on the net [3]. First we converted the original color images (29 images) to gray and then smoothed them by Gaussian kernels by STDs of 1–9.5 with steps of 0.5 (using function cvSmooth of OpenCV). So we have 29 stacks, for which every individual stack consists of 18 images of the same content but with different amounts of blur. Furthermore, every image was distorted with white Gaussian noise by STDs of 5, 10 and 15. So in total we have 29  18  4¼2088 test images (including noise free images). We applied our method, Ducottet’s method and 16 existing methods to this database. For 15 methods, we used the code from [25]. For Ducottet and RHLWT methods, we re-implemented the codes. As said before, if a blur metric is increasing with respect to the blurriness, the calculated amount of blur for every image should be less than the calculated amount of blur of all other images with more real blur in the same

Table 3 Parameters for our method. Noise STD

Smin

Ds

Ns

Nf

smax

0 5 10 15

1 2 3 4

3 3 3 3

4 4 4 4

4 4 4 4

12 12 12 12

Table 4 MAPE for Ducottet’s method and our method. Parameters

SpaDu

Ds Ns Nf Smin ¼ 1 Smin ¼ 2 Smin ¼ 3 Smin ¼ 4

k Noise STD k 0 5 10 15

FftDu

1 4 4

1 4 4

8.2 8.6 11.9 15.6

5.8 6.9 8.2 13.3

SpaDu

FftDu

SpaDu

FftDu

SpaDu

FftDu

SpaQF

FftQF

1 5 5

1 5 5

1 6 6

1 6 6

1 8 8

1 8 8

3 4 4

3 4 4

5.3 7.0 6.1 8.5

4.2 4.8 3.7 8.0

4.3 5.8 3.9 5.8

2.5 3.5 3.8 6.8

3.3 3.2 3.1 5.7

1.8 1.9 4.3 7.9

6.4 6.3 5.7 6.6

1.1 2.1 3.8 5.7

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

stack. Any deviation is considered as an error. The number of errors for every method is shown in Table 5. The parameters of our method were similar to the synthetic images and the Ducottet method was calculated twice: once

1997

with Ns ¼ 6 and the other time with Ns ¼ 8. Other parameters were similar to the synthetic images. Every column of Table 5 is sorted based on the number of errors and the results of our method are shown in bold.

 For noise free images, our method in the spatial   

Fig. 5. Results of our method for synthetic images by blur levels from 1 to 10. Every point represents one test image. In every image all edges have the same blur levels: (a) FFT and (b) spatial.

domain has no error and the ZHANG method has one error. For category with noise of STD of 5, the ZHANG method has 2 errors and ours has 4 errors. For noise of STD 10, the ZHANG method has 14 errors and ours has 19 errors. For noise STD of 15, the spectrum method is the best one with 53 errors and our method has 54 errors, while the ZHANG method has 60 errors.

So our method is in the best or second position in any case. In comparison with Ducottet’s method it is clear that our method has a less number of errors. We also report the performance of every method using Spearman rankorder correlation coefficient which is a non-parametric criterion for monotonicity prediction. The results are shown in Table 6 which are consistent with the number of errors criterion. However, using the number of errors, the locations of errors can be explored as well. The details of error for our method, the ZHANG method and the spectrum method are shown in Tables 7 and 8, respectively. In Table 7 the number of error per stack is shown. For noise STD of 0 and 5, the errors of the ZHANG method just occur in paintedhouse stack but for noise STD of 10 and 15 the errors spread in other stacks. The spectrum method has errors in almost all stacks and for all levels of noise. Our method has no error for noise free images and for noise STD of 5, there is one error in carnivaldolls and dancers stacks and there are two errors for ocean stack. For noise STD of 10, the error occurs in

Table 5 The number of errors for Gaussian blurred images. Noise STD¼ 0 SpaQF Zhang SpaDu_8 FftQF RHLWT SpaDu_6 FftDu_8 FftDu_6 auto lap RIEM grad spectrum wavelet JNBM NRJPEG Marziliano var NSS hp Marichal CPBD

Noise STD ¼ 5 0 1 2 2 3 5 6 8 10 27 31 40 54 56 114 161 224 237 712 1127 1212 1895

Zhang SpaQF FftQF SpaDu_8 FftDu_8 RIEM SpaDu_6 spectrum wavelet FftDu_6 NSS grad auto var RHLWT lap hp Marichal JNBM Marziliano NRJPEG CPBD

Noise STD ¼ 10 2 4 11 28 35 41 46 51 51 74 89 111 111 233 390 528 624 1575 2381 3173 3560 3895

Zhang SpaQF FftQF wavelet spectrum SpaDu_8 FftDu_8 SpaDu_6 RIEM NSS FftDu_6 var auto grad RHLWT hp lap NRJPEG JNBM Marziliano CPBD Marichal

Noise STD ¼15 14 19 31 45 51 53 61 77 107 113 114 235 253 264 412 643 826 2824 3104 3716 3987 4466

spectrum SpaQF Zhang wavelet FftDu_8 FftQF SpaDu_8 SpaDu_6 FftDu_6 NSS RIEM var auto grad RHLWT hp lap NRJPEG JNBM Marziliano CPBD Marichal

53 54 60 65 67 73 83 92 103 147 211 236 348 386 462 557 1013 2263 3500 3535 3685 4466

1998

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

Table 6 Spearman rank correlation coefficients for Gaussian blurred images. Noise STD ¼0 SpaQF Zhang RHLWT SpaDu_8 FftQF SpaDu_6 FftDu_8 FftDu_6 lap RIEM grad spectrum auto JNBM wavelet var NRJPEG NSS Marichal Marziliano CPBD hp

Noise STD ¼5 1.00000 0.99993 0.99986 0.99986 0.99986 0.99957 0.99957 0.99936 0.99794 0.99779 0.99715 0.99601 0.99594 0.94961 0.94014 0.93908 0.91965 0.91161 0.90764 0.82634 0.43616 0.26800

Noise STD¼ 10

Zhang SpaQF FftQF SpaDu_8 spectrum FftDu_8 SpaDu_6 RIEM CPBD NSS FftDu_6 auto Marziliano grad var wavelet NRJPEG JNBM RHLWT lap hp Marichal

0.99986 0.99979 0.99929 0.99658 0.99616 0.99488 0.99388 0.99374 0.99360 0.99153 0.98763 0.97231 0.97153 0.95260 0.94249 0.94029 0.91787 0.88022 0.84171 0.42785 0.37347 0.16135

Zhang SpaQF FftQF spectrum SpaDu_8 FftDu_8 SpaDu_6 Marziliano NSS JNBM CPBD FftDu_6 RIEM var wavelet auto grad NRJPEG RHLWT Marichal hp lap

Noise STD ¼ 15 0.99900 0.99900 0.99801 0.99609 0.99395 0.99388 0.98932 0.98484 0.98086 0.98043 0.97893 0.97787 0.96199 0.94335 0.93993 0.88762 0.84627 0.77688 0.61112 0.50000 0.32714 0.00132

SpaQF spectrum Zhang FftQF FftDu_8 SpaDu_8 SpaDu_6 FftDu_6 JNBM NSS Marziliano wavelet var CPBD RIEM auto grad NRJPEG Marichal RHLWT hp lap

0.99580 0.99523 0.99409 0.99374 0.99345 0.99153 0.98719 0.98342 0.96762 0.96648 0.94385 0.94328 0.94200 0.93374 0.91844 0.74841 0.63204 0.56742 0.50000 0.46771 0.41796 0.34252

Table 7 Details of errors for some better methods. Stacks

Bikes Building2 Buildings Caps Carnivaldolls Cemetry Churchandcapitol Coinsinfountain Dancers Flowersonih35 House Lighthouse Lighthouse2 Manfishing Monarch Ocean Paintedhouse Parrots Plane Rapids Sailing1 Sailing2 Sailing3 Sailing4 Statue Stream Studentsculpture Woman Womanhat Sum

Noise STD ¼0

Noise STD ¼5

Noise STD ¼10

Noise STD ¼15

spectrum

Zhang

SpaQF

spectrum

Zhang

SpaQF

spectrum

Zhang

SpaQF

spectrum

Zhang

SpaQF

1 1 1 2 2 2 2 2 2 2 2 2 2 2 1 2 3 3 2 2 2 1 2 1 2 2 2 2 2 54

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 1 1 2 2 1 2 2 2 2 3 2 2 2 1 2 3 3 1 1 2 1 2 1 2 1 2 2 2 51

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 2

0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 4

1 1 1 2 0 1 2 1 2 2 3 2 2 2 2 3 3 3 3 2 1 2 1 1 2 1 2 2 1 51

0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 2 3 0 2 0 4 0 0 0 0 0 0 0 0 14

0 0 0 0 1 0 0 1 1 0 0 1 0 0 1 1 3 1 2 0 1 0 1 1 1 0 2 0 1 19

2 1 2 3 0 2 1 2 2 0 3 2 1 1 1 3 3 3 4 2 0 1 2 1 2 2 2 3 2 53

0 0 0 2 3 2 2 0 0 4 4 2 2 2 2 5 4 1 3 1 8 2 2 3 0 1 2 1 2 60

2 1 2 0 1 1 0 0 2 4 0 2 3 0 2 2 1 6 2 3 1 2 3 5 3 2 1 1 2 54

more stacks. For noise STD of 15, all the three methods have errors almost everywhere. In Table 8, we show the number of errors for every level of blurriness. As can be seen, the ZHANG method has

more errors for higher blur STDs but our method has more errors at low blur levels. For bikes stack, we show the processing time for our method and Ducottet’s in Table 9. The used CPU for

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

1999

Table 8 Location of errors for some better methods. Blur STD

1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5

Noise STD ¼ 0

Noise STD ¼5

Noise STD ¼10

Noise STD ¼ 15

spectrum

Zhang

SpaQF

spectrum

Zhang

SpaQF

spectrum

Zhang

SpaQF

spectrum

Zhang

SpaQF

0 0 0 0 0 0 0 0 0 0 0 0 31 0 0 0 23 0

0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 31 0 0 0 19 1

0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0

1 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 1 31 0 0 1 17 1

0 0 0 0 0 0 0 0 0 0 0 6 0 2 0 1 0 5

7 4 1 1 0 1 0 0 0 0 0 0 0 0 1 1 1 2

0 0 0 0 0 0 0 0 1 1 0 1 29 0 0 1 19 1

0 0 0 0 0 0 1 0 2 0 0 23 0 8 0 7 0 19

8 6 1 2 2 5 1 1 0 2 2 7 0 2 0 4 4 7

Table 9 Processing time per image for bikes stack. Noise STD

SpaDu (6 Scales) (s)

FftDu (6 Scales) (s)

SpaDu (8 Scales) (s)

FftDu (8 Scales) (s)

SpaQF (s)

FftQF (s)

0 5 10 15

120.30 101.36 72.37 52.39

126.72 104.20 75.49 55.82

234.25 203.73 151.09 112.22

243.37 205.59 153.09 115.54

0.61 0.61 0.60 0.60

9.33 9.00 8.63 8.37

Table 10 Ratio of processing times. k Noise STD k

SpaDu ð6 ScalesÞ SpaQF

FftDu ð6ScalesÞ FftQF

SpaDu ð8 ScalesÞ SpaQF

FftDu ð8 ScalesÞ FftQF

0 5 10 15

197 166 121 87

14 12 9 7

384 334 252 187

26 23 18 14

running programs was Core2 Q9550/2.83 GHz (codes were single-thread) with 8 GB RAM. From Table 9, it is clear that the processing time for our method in spatial domain is about 0.6 s per image which is between 87 and 384 times faster in comparison with the Ducottet method in spatial domain and with Ns ¼ 6 and Ns ¼ 8 (the ratios are shown in Table 10). Ducottet in his paper used a scale step of 0.5. With this value, the method would be even slower than that indicated by Table 9. In Fig. 6 we showed two blurred images from Bikes stack with blur levels of STDs of 1 and 10 pixels along with their segmented results using Ducottet and our method. Both multi-scale methods produce fragmented edges. However, the aim of these methods is to give blur level for every detected edge. 6.3. Motion blurred and JPEG2000 compressed images In this sub-section we report the results of applying the proposed method and 16 other methods to the motion

blurred images and also JPEG2000 compressed images. The 29 reference images of Live database are distorted by motion blur with lengths of 5, 10, 15, 20 and 25 pixels. Then all images (including reference images) were distorted by Gaussian noise with standard deviations of 5, 10 and 15. The used JPEG2000 images are reference images of Live database and their compressed ones with bit rates of 0.030227–3.1539. The Spearman rank correlation coefficients for motion blurred and JPEG2000 compressed images are shown in Tables 11 and 12, respectively. In these tables, for every level of noise, the methods are sorted based on the correlation coefficient. The correlation coefficient of our method (SpaQF in Table 11) for motion blurred images without noise is 0.98 which is a good result. However, some other methods have better performance for noise free images. For noisy images, the performance of our method improves in comparison with other methods. In both cases of noise STD 10 and 15, the correlation

2000

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

Fig. 6. Two images of Bikes stack and their segmented results using Ducottet and our method: (a) bikes image with blur of STD of 1, (b) bikes image with blur of STD of 10, (c) image in (a) segmented by Ducottet (spatial) method, (d) image in (b) segmented by Ducottet (spatial) method, (e) image in (a) segmented by our (spatial) method, and (f) image in (b) segmented by our (spatial) method.

Table 11 Spearman rank correlation coefficients for motion blurred images. Noise STD ¼0 JNBM var spectrum lap grad auto Zhang RIEM RHLWT Marichal SpaQF wavelet CPBD NSS NRJPEG hp Marziliano

Noise STD ¼ 5 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.990 0.980 0.890 0.870 0.740 0.730 0.670 0.550

NRJPEG RHLWT lap var spectrum grad auto Zhang RIEM SpaQF wavelet NSS CPBD Marichal hp JNBM Marziliano

Noise STD¼ 10 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.490 0.440 0.300 0.250 0.240 0.110 0.080 0.070

SpaQF lap var spectrum grad auto Zhang RIEM NRJPEG RHLWT NSS wavelet CPBD Marziliano JNBM Marichal hp

Noise STD ¼15 0.330 0.330 0.330 0.330 0.330 0.330 0.330 0.330 0.320 0.310 0.290 0.290 0.280 0.210 0.190 0.170 0.150

SpaQF lap var spectrum grad auto Zhang RIEM NRJPEG RHLWT NSS wavelet CPBD Marziliano JNBM Marichal hp

0.330 0.330 0.330 0.330 0.330 0.330 0.330 0.330 0.320 0.310 0.290 0.290 0.280 0.210 0.190 0.170 0.150

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

2001

Table 12 Spearman rank correlation coefficients for JPEG2000 images. Noise STD¼ 0 grad Marziliano Zhang RIEM CPBD lap auto var spectrum RHLWT NSS Marichal JNBM SpaQF hp NRJPEG wavelet

Noise STD ¼5 1.000 0.998 0.997 0.997 0.997 0.996 0.995 0.994 0.992 0.991 0.981 0.980 0.973 0.777 0.671 0.487 0.397

grad Zhang RIEM auto RHLWT var lap spectrum Marziliano Marichal NRJPEG JNBM CPBD SpaQF NSS hp wavelet

Noise STD ¼ 10 0.999 0.999 0.998 0.997 0.994 0.994 0.994 0.989 0.985 0.983 0.982 0.980 0.963 0.819 0.809 0.663 0.391

coefficients of our method is equal to the maximum value which is 0.33. The low performance of all methods by increasing the noise suggests that all methods are very sensitive to noise in motion blur estimation. For JPEG2000 compressed images (Table 12), the correlation coefficient of our method for noise free images is 0.777. Although this is an acceptable performance but 13 other methods work better for these images. It can be explained by this fact that, in addition to blurring, the compression produces other types of artifacts like ringing. For JPEG2000 images, when the noise level increases, the performance of our method improves in comparison with others. 7. Conclusions The proposed multi-scale method outperforms the other multi-scale method (Ducottet) in speed and also the accuracy is better. Depending on details of the image, the proposed method is several to tens of times faster. So our method is a promising alternative to use the benefits of multi-scale processing. Other advantages of the new method are first it responds to every line or transition only once at the centroid of the feature and second is that it gives local blur estimation. The third advantage of our method is that it is intrinsically robust against noise. We tested the method using synthetic and real images with three types of blurring in the presence of different levels of noise. The results are promising for Gaussian and motion blurred images. For Gaussian blurring, our method is either in the first or second position in comparison with 16 other methods and for motion blur, it has acceptable results. Because this method gives local blur levels for some detected features it is interesting to check its possibility for image fusion.

Acknowledgement This investigation was supported and funded partially by the FWO Research (FWO Grant No. 3G066908).

lap grad RIEM spectrum auto Zhang var NRJPEG RHLWT JNBM Marziliano CPBD SpaQF NSS Marichal hp wavelet

Noise STD ¼15 0.998 0.998 0.997 0.994 0.994 0.994 0.989 0.984 0.982 0.968 0.949 0.840 0.811 0.646 0.514 0.501 0.372

RIEM auto grad Zhang lap var RHLWT spectrum NRJPEG JNBM Marziliano SpaQF CPBD hp Marichal NSS wavelet

0.998 0.995 0.993 0.991 0.987 0.985 0.983 0.980 0.968 0.962 0.826 0.756 0.713 0.566 0.500 0.462 0.351

References [1] S. Soleimani, F. Rooms, W. Philips, L. Tessens, Image fusion using blur estimation, in: ICIP, 2010, pp. 4397–4400. [2] L. Firestone, K. Cook, N. Talsania, K. Preston, Comparison of autofocus methods for automated microscopy, Cytometry 12 (1991) 195–206. [3] H.R. Sheikh, Z. Wang, L. Cormack, A.C. Bovik, Live Image Quality Assessment Database Release 2, /http://live.ece.utexas.edu/ research/qualityS. [4] R. Ferzli, L.J. Karam, A no-reference objective image sharpness metric based on the notion of just noticeable blur (JNB), IEEE Transactions on Image Processing 18 (4) (2009) 717–728. [5] S. Mallat, S. Zhong, Characterization of signals from multiscale edges, IEEE Transactions on Pattern Analysis and Machine Intelligence 14 (7) (1992) 710–732. [6] C. Ducottet, T. Fournel, C. Barat, Scale-adaptive detection and local characterization of edges based on wavelet transform, Signal Processing 84 (11) (2004) 2115–2137. [7] H.R. Sheikh, A.C. Bovik, L. Cormack, No-reference quality assessment using natural scene statistics: JPEG2000, IEEE Transactions on Image Processing 14 (11) (2005) 1918–1927. [8] Z. Wang, H.R. Sheikh, A.C. Bovik, No-reference perceptual quality assessment of JPEG compressed images, in: IEEE International Conference on Image Processing, vol. 1, 2002, pp. I477–I480. [9] S. Erasmus, K. Smith, An automatic focusing and astigmatism correction system for the SEM and CTEM, Microscopy 127 (1982) 185–199. [10] C.F. Batten, Autofocusing and Astigmatism Correction in the Scanning Electron Microscope, Master’s Thesis, University Cambridge, Cambridge, UK, 2000. [11] X. Marichal, W. Ma, H.J. Zhang, Blur determination in the compressed domain using DCT information, in: IEEE International Conference on Image Processing, vol. 2, 1999, pp. 386–390. [12] N. Zhang, A. Vladar, M. Postek, B. Larrabee, A kurtosis-based statistical measure for two-dimensional processes and its application to image sharpness, in: Proceedings of the Section on Physical and Engineering Sciences of American Statistical Society, 2003, pp. 4730–4736. [13] P. Marziliano, F. Dufaux, S. Winkler, T. Ebrahimi, Perceptual blur and ringing metrics: application to JPEG2000, Signal Processing: Image Communication 19 (2) (2004) 163–172. [14] D. Shaked, I. Tastl, Sharpness measure: towards automatic image enhancement, in: Proceedings of IEEE International Conference on Image Processing, vol. 1, 2005, pp. 937–940. [15] R. Ferzli, L.J. Karam, J. Caviedes, A robust image sharpness metric based on kurtosis measurement of wavelet coefficients, in: Proceedings of the First International Workshop on Video Processing and Quality Metrics for Consumer Electronics, 2005. [16] R. Ferzli, L.J. Karam, A no reference objective sharpness metric using Riemannian tensor, in: The Third International Workshop on Video Processing and Quality Metrics for Consumer Electronics VPQM-07, Scottsdale, Arizona, 2007.

2002

S. Soleimani et al. / Signal Processing 93 (2013) 1988–2002

[17] N. Narvekar, L.J. Karam, A no-reference perceptual image sharpness metric based on the cumulative probability of blur detection, First International Workshop on Quality of Multimedia Experience (QoMEX), 2009. [18] J. Kautsky, J. Flusser, B. Zitov, S. Simberov, A new wavelet-based measure of image focus, Pattern Recognition Letters 23 (14) (2002) 1785–1794. [19] P. Perona, J. Malik, Detecting and localizing edges composed of steps, in: Proceedings of the Third International Conference on Computer Vision, 1991, pp. 52–57. [20] P. Kube, P. Perona, Scale-space properties of quadratic feature detectors, IEEE Transactions on Pattern Analysis and Machine Intelligence 18 (10) (1996) 987–999. [21] T.C. Folsom, R.B. Pinter, Primitive features by steering, quadrature, and scale, IEEE Transactions on Pattern Analysis and Machine Intelligence 20 (11) (1998) 1161–1173.

[22] W.T. Freeman, E.H. Adelson, The design and use of steerable filters, IEEE Transactions on Pattern Analysis and Machine Intelligence 13 (9) (1991) 891–906. [23] M. Frigo, S.G. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93 (2) (2005) 216–231. (special issue on ‘‘Program Generation, Optimization, and Platform Adaptation’’). [24] J.M. Geusebroek, A.W.M. Smeulders, J. van de Weijer, Fast anisotropic Gauss filtering, IEEE Transactions on Image Processing 12 (8) (2003) 938–943. [25] A.V. Murthy, L.J. Karam, Ivquest and Video Quality Evaluation Software, /http://ivulab.asu.edu/Quality/IVQUESTS.