Analysis of the Image Spectrum for Distortion Diagnostics

The diagnostic of possible distortions in given visual data is a significant problem in signal and image processing. In the paper, we discuss the task of automatic determination of the type and parameters of expected distortion on the analysis of the obtained image. In solving the problem, we restrict the research by the identification of linear homogeneous convolution operators of the following three types: circular form with rectangular profile, circular form with Gaussian profile, and linear motion blur ones. The diagnostics identification procedures base on the image spectrum analysis and use the mean radial profile of the image spectrum amplitude. For every distortion type, we elaborate corresponding diagnostics procedures, that when combined into one common process, form a unique recognition algorithm. The studies and experiments demonstrate the reliability of distortion identification algorithms against signal precision and noise influence.


Introduction
The quality of the registered image depends on the overall optical system parameters and the shooting process. As a rule, the distortions of a signal are the results of focusing inaccuracy, relative movement of a camera and observed scene, radiation scattering in transit medium, non-ideal aperture function. If the camera characteristics, position, movement, optical properties of media, etc. are known, we can calculate the point spread function and decide the image restoration task. Unfortunately, in practice, all these data are usually unknown. In many cases, we may only suppose the class of possible image defects, while the real type and distortion parameters, that are necessary to form the reconstruction operator, remain unknown. These data can be established by studying the obtained signal only. Therefore, the diagnostic of distortion type and parameters by means of observed image analysis is important and actual.
There are many studies concerned the decision of this problem. Their main part dedicates to the estimation of linear smoothing parameters, appearing due to the relative shift of the observable scene and registration device [1]- [3]. Significantly lesser part orients to circular blur parameters determination [4]- [7], including anisotropic one [8], and a small part of investigations intend for Gaussian blur estimation [9]. In spite of plenty of publications concerning these questions, the problem of distortion type and parameters diagnostic is insufficiently studied.
We restrict this paper to the questions for diagnostics the type and parameters of the homogeneous linear distorting operator. We suppose that detectable distortions comply with one the following twodimensional functions: circular one with a rectangular profile (rectangular blurring), circular one with Gaussian profile (Gaussian blurring), and linear one with a rectangular profile (linear smoothing). First Here g(u,v) is an ideal signal of the observed scene, h(x,y,u,v) is distorting operator, γ{·} defines system nonlinearities, and ξ is additive non-correlated noise. If impulse response is space-invariant, then h(x,y,u,v) is the function of differences of arguments: h(x-u,y-v). If also system nonlinearities γ{·} may be neglected, we obtain the convolution integral: Here h(x,y) is the spatially invariant point spread function or, in other words, the smoothing aperture. According to the convolution theorem, equation (1) may be written as the product of spectrums of signal G(ω) (ω = {ω x ,ω y }) and of scattering function H(ω) in the combination of noise spectrum Ξ(ω): (2) The amplitude of the real image spectrum does not possess any specific properties [10]. Hence, only very general conclusions about the image characteristics may be deduced basing on the form of image spectrum G(ω) [11], [12]. The spectrum Ξ(ω) of independent uncorrelated noise may be assumed as uniform. From the other side, distorting functions spectra H(ω) are unique and noticeably differ from real images spectra. According to (2), the distinguishing features of H(ω) are transferred to corrupted image spectra F(ω) [13], [14]. So, if H(ω) possess some specific amplitude-frequency characteristics, then these features may be detected in F(ω) spectrum and used for distorting operator diagnostics.
The distorting operator is defined by the kernel h(x,y) in (1). In real tasks, h(x,y) is the function with finite support, hence h(x,y) = 0 at (|x|,|y|) ≥ d h . From the other side, an image of a real scene extends to the whole image plane and is spatially unlimited function. It follows that if the distorting function h(x,y) is spatially limited, then the spectrum of the resulting image will have ordered arrangement of zeroes. Therefore, once distinctive features of F(ω) are found, one might recognize the type and the parameters of distorting operator H(ω).

Radial profile of spectrum amplitude
Let be restricted by diagnostics of the linear space-invariant distortions described by equation (1). Two of three concerned below distorting operators have circular symmetry. Hence, every central section of the two-dimensional signal spectrum holds identical information content concerning the distortion. Along with previously mentioned, every section also contains some part of the image signal information associated with the direction of the section. This information should be suppressed. To analyze the average alteration of spectrum amplitude with the frequency, form the following function.
Let M(ω,φ) be the spectrum amplitude of F(ω) in polar coordinates (ω,φ). Averaging the values of M(ω,φ) over φ we obtain the mean spectrum amplitude value for the frequency ω: In the discrete case of this transformation, P(ω) may be calculated as (ω) (ω) (ω ,ω ) / (ω) Here C(ω) is zero-centered ring of radius ω and width Δω containing spectral points (ω x ,ω y ) χ C(ω); M(ω x ,ω y ) is the spectrum amplitude at the point (ω x ,ω y ), and N(ω) is the number of pixels hitting the ring C(ω). Let name the function P(ω) as the mean radial profile of image spectrum amplitude. As experiments demonstrate, P(ω) is rapidly decreasing with increasing ω, and the form of P(ω) enough little depends on the content of "good quality" image.

Circular blurring
The distorting operator of circular blurring with a rectangular profile and radius r is the following: Two-dimensional operator h C (x,y) in its section is rectangular smoothing impulse. The amplitude characteristic of the spectrum of one-dimensional square-wave and unit in ±r impulse is: S(ω) = 2r |sin(ωr)|/(ωr). (6) This function equals zero in the points ω = ±ρn, where ρ = π/r, n = 1, 2, ... . These nulls will be present in the two-dimensional spectrum too, forming concentric rings of corresponding radii. Therefore, the determination of distortion consists in finding the radius period ρ of these null-rings and calculating the scattering radius r. The following cepstrum transformation [15], [16] may assist in this analysis: Here |·| means signal amplitude and is inverse Fourier transform. Cepstrum (7) may be treated as the power spectrum of the function log |S(ω)| 2 . Figure  The distortion diagnostics may be realized as with Fourier spectrum magnitude, shown in figure 1(c), as with the result of cepstrum transformation (7). The experiments demonstrate that the second variant is significantly more sensitive to data inaccuracy and noise presence. Spectrum magnitude and it's mean radial profile (4) are more stable to distortion diagnostics.
The ring structure in the smoothed image spectrum may be analyzed with the mean radial profile of image spectrum amplitude (4). The period ρ of minima/maxima, existing in profile P(ω), may be found as the point of maximum of the following functional: Here W is the image size, and [·] means the integer part of the argument. If the variation of P(ω) is less than the range of P(ω) plus some small threshold, then the algorithm concludes that circular blurring does not detect. Otherwise, when ρ C is found, the blurring radius r C will be:  figure 1(a); b) smoothed image in figure 1(b); c) precision reduction and Gaussian noise addition: 0noiseless, 3, 5, and 10 -noise σ = 3, 5, and 10 gradations.

Gaussian blurring
The zero-centered distorting operator of Gaussian blurring with circular symmetry is the following: h G (x,y) = C exp{-(x 2 + y 2 )/2σ 2 }. (10) Function h G (x,y) is not spatially limited and in any of its section has the profile of normal distribution. According to (2), the Fourier spectrum of an image distorted by Gaussian blurs will be: F(ω x ,ω y ) = C exp{-(x 2 + y 2 )/2σ 2 }G(ω x ,ω y ) + Ξ(ω x ,ω y ). Therefore the mean radial profile of image spectrum amplitude (3) for the distorted image is: Here G M (ω) and Ξ M (ω) are the mean values of spectrum amplitude for spatial frequency ω of observable scene spectrum and additive uncorrelated random noise.
With increasing ω, the first item in (11) decreases enough fast, while the second one may be considered as invariable. In reality, the power of noise does not exceed several percents of signal power. At that, the ω-space characterizes as containing two areas of different ratios of summands. For amplitude P(ω) at the range of large ω, when C exp{-ω 2 /2σ 2 }G M (ω) < Ξ M (ω)/2, we may write: log P(ω) ≈ log Ξ M (ω) ≈ C 1 , That is, the log P(ω) values at large ω are close to a constant definable by the noise level Ξ M . On the other side, in the range of small ω, i.e. when C exp{-ω 2 /2σ 2 }G M (ω) > 2Ξ M (ω), log P(ω) ≈ -ω 2 /2σ 2 + log G M (ω) + C 2 = A(ω) + B(ω) + C. (12) Here the first component A(ω) = -ω 2 /2σ 2 is defined by Gaussian blur, while the second one B(ω) = log G M (ω) depends on the specified image, and C is some constant.
Gaussian blur estimation reduces to finding the parameter σ which defines A(ω) component in (12); that is we need to determine the parameter of someone signal from the mixture of two. The component B(ω) depends on the specific image. The analysis of real images demonstrates that in the relevant range of ω, the component B(ω) is convex downwards with the exception of a tiny area of ω values near zero (see figure 2(a) ). From the other side, A(ω) is described by the formula y = -ax 2 + b and is convex upwards everywhere. The experiments with the images distorted by the operator (10) reveal that in the considered section of ω, the curvature index of A(ω) is greater than one of B(ω) component (in absolute values). As a result, it is possible to find the range (ω x ,ω y ) in summary curve (12), where log P(ω) should be convex upwards.
This means that the plot of the mean radial profile of image spectrum amplitude log P(ω), calculated for an image distorted by Gaussian circular blur (10), should contain some range, which is close to parabola y = -ax 2 + b with a = 1/(2σ 2 ). The values a and b can be obtained with the help of Hough transform, consisting in mapping a set of analyzed points to the space of estimated function parameters, and generating the accumulation diagram which corresponds to supposed analytical dependence [11], [17]. Therefore, if the plot of the mean radial profile contains a region appropriate to parabola -ax 2 + b, then noticeable nodal point (extremum of the distribution) will appear in Hough transformation diagram with coordinates corresponding to required values a and b. If a value is found, the dispersion value of Gaussian blur will be calculated as follows: (13) The choice of the range of mean radial amplitude profile, to be used in constructing the Hough transformation diagram, is the important stage. For this purpose, we draw the tangent line from below to initial part of the obtained curve. The tangent line is illustrated by the biased line 'T' in figure 3(b). The contact points ω 1 and ω 2 restrict the range for analysis. If only one point is detected (i.e. log P(ω) is convex downwards), or ω 1 and ω 2 are very close, one may decide that Gaussian circular blurring is not detected.
Note that a value given by the proposed method will be distorted by non-eliminated component B(ω) in (12). Its action will reduce the log P(ω) curvature and hence will result in the increased a value. In this connection, it seems perspective to estimate the average shape of log G M (ω) component on the base of the representative set of good quality images, and then to subtract the obtained averaged component from the right part of eq. (12) for the range (ω 1 ,ω 2 ). a) b) c) Figure 3. Plots of the mean radial profile of spectrum after Gaussian blurring (10) with σ=5 pixels: a) full image precision; b) image precision reduced to 256 gradations (line 0) and Gaussian noise of 1 and 2 gradations RMSE is added (lines 1 and 2). c) The accumulation diagram in Hough transformation space for the region (ω 1 ,ω 2 ) of profile 0 in (b).

Linear blurring
Linear blurring is the result of linear monotonic relative displacement between the observed scene and registering device. Its parameters are the direction α (0 ≤ α < 2π) and the extension r of the displacement. The centrally symmetrical linear blurring operator h(x,y) is described by the formula: This operator specifies the convolution of an image signal with the segment of length r directed at the angle α. The spectrum amplitude of the operator (14) is S(ω x ,ω y ) = 2|sin(Ωr/2)|/(Ωr/2), where Ω = 2 2 ω ω cos(α arctan(ω ω )) Therefore, S(ω x ,ω y ) looks like a set of strips oriented orthogonally to blurring direction α and are positioned with the period W/r, where W is the image size and r is the extension of the blur. The section profile of S(ω x ,ω y ) in the direction α will be |sin x|/x. The Radon transform is usually used for diagnosis of such distortions [2]. We propose to use the values C(q 1 ,q 2 ) calculated with cepstrum transformation (7). For an image distorted by linear blur, the periodic structure (strips of zeroes) will appear in its spectrum, and its cepstrum will contain zero secant line of increased values with a set of local maxima. This line will be orthogonal to spectrum strips and will coincide with blurring direction. Such a signal may be detected by the analysis of the cepstrum in polar coordinates (α,q), where α is the angle and q is the center distance. The required angle α Res corresponds to the α-directional maximum of the integral along zero transient line C S (α,q): If the value of the integral in the right part of the equation (16) in direction α Res does not exceed the α-mean value of this integral by some given threshold, then we consider that linear blurring is not detected.
The local maxima along the found secant line (α Res ,q), mentioned above, arrange with the period r equal to blurring size. Theoretically, this fact makes possible the simultaneous determination of the direction and extension of the blur. However, due to appreciable fluency of noise, these extrema will be imperceptible; therefore, they may not serve as a reliable base for extension estimation. Below we discuss a substantially more stable method for strips period determination, which, per se, is equivalent to the analysis of single respective projection from Radon transformation [18].
Let (α Res ,q) be zero-transit straight line oriented in blur direction α Res . Through every point q in (α Res ,q) let draw the line L ψ,q (l) at the angle ψ = α Res + π/2, orthogonal to α Res . For every q calculate the average value P α (q) of the spectrum amplitude logarithm log |S(ω x ,ω y )| along this line (ω x ,ω y ) χ L ψ,q (l) and obtain the mean profile for α Res direction: Here L means the length of L ψ,q within the borders of spectrum space; ω x and ω y are calculated as follows: ω x = q cos α + l sin α and ω y = q sin α -l cos α. Then, analogously to (8), on already known P α (q) find the strip period ρ L as the point giving the maximum of the functional: When the period value ρ L and the image size W are known, the blurring extension value r S founds as:  figure 1(a) by linear blur of 32 pixels and 60° rotation; b) spectrum of (a) in logarithmic scale; c) mean profile P α (q) of spectrum (b) amplitude; d) cepstrum of (a); e) spectrum of the image in (a), distorted by Gaussian noise σ=7 gradations; f) mean profile P α (q) of spectrum (e) amplitude.  Figure 4 illustrates the analysis process; figure (a) is the image in figure 1(a) distorted by blurequivalent signal (blur extension 32 pixels, rotation 60°). The spectrum (in logarithm scale) of the image in figure 4(a) in precise form is in (b), and of the image in the reduced precision form (256 scales) with added Gaussian noise of σ = 7 scales RMSE is in (e). The figure (d) demonstrates the cepstrum (7) of (b). Plots in figures (c) and (f) represent the dependences P α (q) in (17) and they are the mean radial profiles of image spectrum amplitude for figures (b) and (e) correspondingly.

Conclusion
In the paper, we studied the distortion assessment problem basing on observed image analysis, when no supplementary information is present nor about the type and parameters, nor even about the presence of any distortion. Proposed algorithms demonstrate their capabilities for estimating the type and the parameters of following linear homogeneous distortion operators: circular blur with rectangular profile, circular blur with Gaussian profile, and linear motion blur. The experiments demonstrate enough high precision of blur parameters identification against the signal quantization and noise addition.