Improving the accuracy of restoring a distorted image via determining the distortion parameters from the Fourier spectrum

Abstract

The spectral method is developed for estimating the parameters of the point-spread function (PSF) in the problem of restorating the distorted (smeared, defocused) images. The method is based on the analysis of the Fourier spectrum of a distorted image. This method makes it possible to estimate the PSF parameters: the angle θ and magnitude Δ of image smearing, as well as the size r of image defocusing spot. New estimates are obtained for parameters θ and Δ using the Nyquist frequency and a estimate for r using the Bessel function. The results of applying this method to image processing are presented. This method can be used to enhance the accuracy of smeared and defocused image restoration via their mathematical processing by stable methods for solving the Fredholm integral equations of the first kind (ill-posed problem).

Keywords

• Image distortions (smearing
• defocusing)
• Point-spread function
• Distortion parameters
• Fourier spectrum
• Integral equations
• MatLab

-

Various image recording devices (IRDs), namely, digital photo cameras, video cameras, tracking systems, telescopes, microscopes, et al. record images of objects – people, animals, technical details, car license plates, astronomical objects, biological microorganisms, etc. In this case, an image may be smeared (due to the device shift or the object movement) or defocused (due to an erroneous focus setting), as well as noisy by external or instrumental noise [ 1–4 ]. These image distortions can be eliminated by a technical way – via photography at a stationary IRD or object (in this case, there will be no image smearing), via correct setting the focus (there will be no defocusing), and via photography in the absence of noise.

However, photography conditions are not always favorable. Examples: photographing a fast-flying aircraft (the photograph may become smeared) [4, p. 172], fuzzy initial pictures of astronomical objects taken by the Hubble telescope [4, p. 105], [ 5 ], etc.

The distorted images can be restored via mathematical and computer processing ([ 2, 4, 6-11 ], etc.). Such a way is an essential addition to the image restoration problem, when technical restoration is problematic. 2

Mathematical formulation of the image restoration problem The problem for eliminating of image smearing or defocusing is usually realized by solving two-dimensional Fredholm integral equation (IE) of the first kind of convolution type [ 2, 4, 7, 9–12 ]: ¥ ¥ ò ò h(x - x, y - h) w(x,h) dx dh = g(x, y) + dg , -¥ -¥ where h is the PSF or the kernel of IE, which is usually spatially invariant (a difference function); w and g are the intensity distribution over the true and distorted images, respectively; δg is the noise.

The goal of the work is to enhance the accuracy of the solution of IE (1) via refining the PSF, or the kernel h(x, y) of IE (1). 3

Methods for solving integral equation (1) The problem for solving IE (1) is ill-posed (essentially unstable) [ 4, 7–11 ]. Therefore, we will use for its solving stable methods – the Tikhonov regularization [ 2, 4, 7–13 ], as well as the Wiener parametric filtering [ 2, 4, 9 ].

The solution of IE (1) by Tikhonov’s regularization with Fourier transform (FT) is wa (x, y) = F -1(Wa (w1, w2)) , where F -1 is the inverse FT, a > 0 is the regularization parameter, and Wa (w1, w2) is the regularized spectrum (two-dimensional FT) of the solution:

Wa (w1, w2 ) =

H *(w1, w2 ) G(w1, w2 ) H (w1, w2 ) 2 + a (w12 + w2 ) p 2 (1) (2) where p ³ 0 is the regularization order (usually p = 1, 2 or 3), spectra (FTs) H (w1, w2) = F(h(x, y)), G(w1, w2) = F(g(x, y)) , and F is the direct FT.

The solution of IE (1) by Wiener’s parametric filtering is wK (x, y) = F -1(WK (w1, w2)) , where K ³ 0 is the parameter equal to the noise/signal power ratio (NSPR) and the solution spectrum equals WK (w1, w2 ) =

H *(Hw(1w,w1,2w)2G)(2w+1,Kw2 ) . (3) The solution of IE (1) by Tikhonov's regularization (TR) and Wiener’s parametric filtering (WPF) is realized in the IPT package of the MatLab system in the form of m-functions deconvreg.m and deconvwnr.m, respectively. However, even the TR and WPF methods are very sensitive to errors of the distortion parameters – the smearing values Δ and θ, as well as ρ of image defocusing, i.e. to inaccuracies in knowledge of the PSF, or to errors of the kernel h of IE (1).

As an illustration, Fig. 1a shows the image of the phantom tomogram developed by us [4, p. 14] mrt-1-02.bmp 407×380 pixels, smeared at angle θ = 35°, smear Δ = 14 pixels.

Fig. 1b shows the image restoration result via solving IE (1) by the TR method (TRM) according to (2) with a = 10–4 and the WPF method (WPFM) according to (3) with K = 10–4 for exact smear parameters θ = 35° and Δ = 14 pixels. Furthermore, the m-functions of the system MatLab fspecial.m, imfilter.m (direct problem), as well as deconvreg.m and deconvwnr.m (inverse problem) are used [ 6 ]. And Fig. 1c shows a restored image with erroneous values of the smear parameters ~q = 37° and D = 16 pixels (a = K = 10–3). The values of the pa~ rameters a and K are chosen by the selection way, namely, via displaying images with various parameters a and K and visual selection of the parameters’ ~ values giving the best image restoration wa and wK . Although q and D~ differ insignificantly from the exact θ and Δ, the restoration is unsatisfactory (Fig. 1с).

Fig. 2a presents a similar example with colored (rgb) defocused image of the astronomical object (galaxy) m83.jpg 378×400×3 pixels.

In this example, the PSF is a uniform disk of radius ρ = 10 pixels. Fig. 2b shows the image restoration result via solving IE (1) by the TR and WPF methods with exact defocusing parameter ρ = 10 pixels (a = K = 10–4). And in Fig. 2c, a restored image is given with erroneous value of the parameter ~r = 11 pixels (a = K = 2·10–3). Although ~r differs little from the exact ρ, the restoration is unsatisfactory (Fig. 2c). This example is given for the comparison with the results of mathematical image restoration obtained by the Hubble telescope [ 5 ] in the presence of remanence spherical aberration of the telescope mirror which is equivalent to defocusing.

These and other examples [ 4, 7, 10, 11, 14, 15 ] indicate that some method is necessary to enhance the accuracy of the estimation of image distortion parameters, in other words, to enhance the accuracy of knowledge of the point-spread function. One of these methods is the spectral method for estimating the PSF [ 4, 10, 11, 14 ] which determines the PSF parameters by Fourier spectrum of the distorted image. In this paper, the spectral method is further developed, in particular, its verification is performed on a number of distorted images.

Note the following existing ways for estimating the image distortion parameters: from streaks in the image [ 4, 9 ], from the blurring of points on the image in the case of defocusing [ 2, 4, 9 ], from the Fourier spectrum of the image [ 1, 16, 17 ] and others. We should point out the rapidly developed methods of “blind” and “semiblind” deconvolution [ 2, 6, 18 ] for estimating the PSF. We should also point out [ 7, 19 ], in which a stable algorithm was developed for reconstructing images under inexactly known response function (or PSF). 4

Spectral method Denote by g(x, y) the intensity distribution over the distorted image, where the x axis is directed horizontally and y vertically down.

We carry out a two-dimensional Fourier transform (FT) of distorted image g(x, y) (4) where wx and w y are Fourier frequencies directed horizontally and vertically, like x and y. We suppose that the FT (4) is calculated through discrete FT (DFT) with the help of m-function fft2.m. Moreover, the DFT centering procedure [6, p. 126] is performed with the help of fftshift.m. As a result, we obtain the complex DFT (the Fourier spectrum) G(wx , wy ) which is conveniently expressed as Re G(wx , wy ) or modulus | G(wx , wy ) |. 4.1

Estimating the image smear parameters Let us consider smeared image of the phantom tomogram. (Fig. 1a). We introduce new axes on the smeared image, namely, we direct axis x' along the smear and axis y' perpendicular to the smear. On the spectrogram, we introduce w axis along the smear. As a result of image smear along x', the Fourier spectrum G(wx , wy ) shrinks along w axis (the decrease of the wmax occurs). This decrease grows with the raise of the smear value Δ. Such effect is connected with suppression of the high Fourier frequencies when image smearing [ 1, 16, 17 ].

Fig. 3 shows three versions of the Fourier spectrum in the form of centered discrete FT (DFT) of the image mrt-1-02.bmp. Fig. 3a represented the spectrum of the undistorted image, Fig. 3b and Fig. 3c shows the spectrum of the smeared image Re G(wx, wy ) and | G(wx , wy ) | respectively. One can see that the spectrum appearance is substantially different. This allows us to determine whether the image was smeared initially. Furthermore, comparison of Fig. 3b and Fig. 3c shows that the spectrum modulus | G(wx , wy ) | gives sharper image than Re G(wx, wy ) . age Re G(wx ,wy ) , c – DFT of the smeared image | G(wx ,wy ) | .

Estimating the image smear magnitude. The smeared image spectrum (Fig. 3c) consists of a set of almost parallel lines and a central quasi-ellipse (сf. [ 1–4, 10, 11, 14, 16, 17 ]). Let us draw through the quasi-ellipse the middle line L and the axis w perpendicular to it, as well as the horizontal axis wx and the vertical axis w y . We also note on the axis w the first (when w = w1 ) and the last (when w = wmax ) zeros of the spectrum | G(w) |. As shown in [ 11 ], smear magnitude equals

D = 2 wwm1ax . (5) In relation to the considered example, the dimensionless ratio wmax w1 (in any same units: pixels, centimeters, Nyquist frequencies, etc.) is estimated from several measurements in Fig. 3c as 7.02 ± 0.08 pixels. Therefore, D = 14.04 ± 0.16 pixels, which is close to the accurate value of D = 14 pixels.

In [ 10, 14 ], to determine Δ, the formula: D = 2 w2 w1 is used (see Fig. 3c). However, formula (5) gives a more accurate result (see also [ 11 ]). Estimating the smear direction. Using Fig. 3c we determine the angle ~q between the horizontal axis wx and w axis, as well as the angle y~ = 90° - ~q (measured angles). However, the angles ~q and y~ , generally speaking, do not coincide with the true angles θ and ψ. This is due to the fact that the image in Fig. 1 and its spectrum in Fig. 3 are, generally speaking, rectangular in size M ´ N (M rows and N columns). Their ratio is r = M N . Only when M = N and ~ therefore r = 1 (when the image is a square) q = q and y~ = y .

To determine the true angles θ and ψ we take into account that tg y equals to the slope coefficient of any straight line, and when stretching or compressing the image, i.e. when r is changed, the slope coefficient changes r times: tg y~ = r × tg y , and therefore the true angle ψ equals ae tg y~ ö y = arctg ç ÷ è r ø (6) and q = 90° - y .

In our example, according to Fig. 3, we determine by several measurements: ~ ~ q » 33°.0 ± 0.4, y~ = 90° - q » 57°.0 ± 0.6, M = 407 , N = 380 , r = M N = 1.071. Using (6), we get: y = 55°.2 ± 0.6 , and q = 34°.8 ± 0.4 , which is close to the exact value of smear angle q = 35° .

So, using the spectral method, we determined with good accuracy the smear magnitude (an integer): D = 14 pixels and smear angle: q = 34°.8 ± 0.4 » 35° which are close to the exact values.

Now, using the values Δ and θ found from the spectrum, we can obtain a high-quality restoration of the tomogram–phantom image via solving IE (1) by Tikhonov's regularization and Wiener’s parametric filtering using m-functions deconvreg.m and deconvwnr.m. In this case, the point-spread function (PSF) was calculated with the help of m-function fspecial.m [ 6 ]. The result of restoration is shown in Fig. 1b. 4.2

Estimating the image defocusing parameter Consider the simplest variant of defocusing, when every point on the object is converted in its image into a uniform circle (disk) of radius ρ and density 1/πρ2 [10, p. 158]. This can occur in the case of a thin lens with a circular aperture [4, p. 193]. Consider one such circle. Two-dimensional Fourier transform of a uniform circle of radius ρ (its optical transfer function – OTF) is expressed through the one-dimensional Hankel transform [20, p. 69], [21, p. 249]: F (w1,w2 ) = òòei(w1x+w2 y) dx dy == p2rp2 rò J0 (wr) r dr = w22r2 wòrJ0 (z) z dz , D 0 0 (7) where D is the circle area, w = w12 + w22 , and J0(z) is the first-kind Bessel function of the zeroth order. The last integral in (7) equals [21, p. 668] wr ò J0 (z) z dz = wr J1(wr) , (8) 0 where J1(ωρ) is the first-kind Bessel function of the first order. Taking into account (8), we obtain (cf. [16, p. 24], [17, p. 100]):

F (w) = w2r J1(wr). (9) r = 3.84 w1, 7.02 w2, 10.16 w3, 13.32 w4 ,!, (11) where w1, w2, w3, w4, ! are the Nyquist frequencies (but not in pixels) corresponding to the semiaxis of each ellipse. Upon discretization, the maximum Nyquist frequency is ωmax = π both along the horizontal and vertical axes in Fig. 5c. Then the frequencies wi , i = 1, 2, 3, 4,… are equal to wi = wi rel wmax = wi rel p, i = 1, 2, 3, 4,!, where wi rel = wi wmax are dimensionless ratios.

We obtain: w1rel = 0.122 , w1 = 0.383 , r = 10.02 ; w2rel = 0.226 , w2 = 0.710 , r = 9.88 , etc. On the average, r = 9.95 ± 0.07 , which is close to the exact value of r = 10 pixels. cused image Re G(wx ,wy ) , c – DFT of the defocused image | G(wx ,wy ) | . Only now, using the value r = 9.95 found from the spectrum and rounded to r = 10 , we got the opportunity to restore with increased accuracy the image of the galaxy M33 via solving IE (1) by the TR (and WPF) method according to (2) with a = K = 10-4 using m-functions deconvreg.m and deconvwnr.m. The restoration result is shown in Fig. 2c. We see a clear image restoration and it is due to the fact that the spectral method makes it possible to determine the defocusing parameter ρ almost exactly.

We note that in papers [ 4, 10, 11 ] the defocusing variant is also considered in the case when the PSF is a Gaussian, and the estimates of the parameter s of the Gaussian are used. 5

Conclusion This work confirms the effectiveness of the spectral method for determining the image distortion type (smear or defocusing) and determining the distortion parameters. The proposed spectral method can be applied to enhance the resolution of image recording devices (digital video cameras, telescopes, microscopes, tomographs, etc.).

This work was supported by the grant MFKTU ITMO (Project No. 619296).

References

1. Bates , R.H.T. , McDonnell , M.J. : Image Restoration and Reconstruction . Oxford U. Press, Oxford ( 1986 ).

2. Gonzalez , R.C. , Woods , R.E. : Digital Image Processing. 2nd edn. Prentice Hall , Upper Saddle River ( 2002 ).

3. Jähne , B. : Digital Image Processing . Springer, Berlin-Heidelberg ( 2005 ).

4. Sizikov , V.S. : Direct and Inverse Problems of Image Restoration, Spectroscopy and Tomography with MatLab . Lan' Publ., St. Petersburg ( 2017 ) (in Russian) .

5. Saïsse , M. , Rousselet , K. , Lazarides , E.: Modeling technique for the Hubble Space Telescope wave-front deformation . Applied Optics 34 ( 13 ), 2278 - 2283 ( 1995 ).

6. Gonsales , R.C. , Woods , R.E. , Eddins , S.L. : Digital Image Processing using MATLAB . Prentice Hall, New Jersey ( 2004 ).

7. Voskoboinikov , Yu .E.: A combined nonlinear contrast image reconstruction algorithm under inexact point-spread function . Optoelectronics, Instrumentation, and Data Processing 43 ( 6 ), 489 - 499 ( 2000 ). doi: 10 .3103/S8756699007060015.

8. Hansen , P.C. : Discrete inverse problems: Insight and algorithms . SIAM, Philadelphia ( 2010 ).

9. Sizikov V.S. , Ékzemplyarov R.A. : Operating sequence when noise is being filtered on distorted images . J. Optical Technology 80 ( 1 ), 28 - 34 ( 2013 ). doi: 10 .1364/JOT.80.000028.

10. Sizikov , V. S. : Spectral method for estimating the point-spread function in the task of eliminating image distortions . J. Optical Technology 84 ( 2 ), 95 - 101 ( 2017 ). doi: 10 .1364/JOT.84.000095.

11. Sizikov , V.S. , Stepanov , A.V. , Mezhenin , A.V. , Burlov , D.I. , Éksemplyarov , R.A. : Determining image-distortion parameters by spectral means when processing pictures of the earth's surface obtained from satellites and aircraft . J. Optical Technology 85 ( 4 ), 203 - 210 ( 2018 ). doi: 10 .1364/JOT.85.000203.

12. Gorshkov , A.V. : Improvement of image resolution during physical experiment data processing and determination of unknown test response by the REIMAGE packet programs . Instrumentation and Experimental Technique 2 , 68 - 78 ( 1995 ).

13. Tikhonov , A.N. , Goncharsky , A.V. , Stepanov , V.V. , Yagola , A.G. : Numerical Methods for the Solution of Ill-Posed Problems . Kluwer, Dordrecht ( 1995 ).

14. Sizikov , V.S.: Estimating the point-spread function from the spectrum of a distorted tomographic image . J. Optical Technology 82 ( 10 ), 655 - 658 ( 2015 ). doi: 10 .1364/JOT.82.000655.

15. Donatelli , M. , Huckle , T. , Mazza , M. , Sesana , D. : Image deblurring by sparsity constraint on the Fourier coefficients . Numerical Algorithms 72 ( 2 ), 341 - 361 ( 2016 ).

16. Vasilenko , G.I. , Taratorin , A.M. : Image Restoration. Radio i Svyaz' Publ., Moscow ( 1986 ) (in Russian) .

17. Gruzman , I.S. , Kirichuk , V.S. , Kosykh , V.P. , Peretyagin , G.I. , Spektor , A.A. : Digital Image Processing in Information Systems . NSTU Publ., Novosibirsk ( 2002 ) (in Russian) .

18. Yan , L. , Lin , H. , Zhong , S. , Fang , H.: Semi-blind spectral deconvolution with adaptive Tikhonov regularization . Applied Spectroscopy 66 ( 11 ), 1334 - 1346 ( 2012 ).

19. Voskoboinikov , Yu .E., Mukhina , I.N. : Local regularizing algorithm for high-contrast image and signal restoration . Optoelectronics, Instrumentation, and Data Processing 43 ( 3 ), 41 - 48 ( 2000 ).

20. Bracewell , R.N.: The Hartley Transform . Oxford University Press, New York ( 1986 ).

21. Korn , G.A. , Korn , T.M. : Mathematical Handbook for Scientists and Engineers . McGraw-Hill , New York ( 1961 ).

22. Donatelli , M. , Estatico , C. , Martinelli , A. , Serra-Capizzano S.: Improved image deblurring with anti-reflective boundary conditions and re-blurring . Inverse Problems 22 , 2035 - 2053 ( 2006 ). doi: 10 .1088/ 0266 -5611/22/6/008.