CN104820969A - Real-time blind image restoration method - Google Patents

Abstract

The invention provides a real-time blind image restoration method, which comprises the following steps of carrying out Fourier transformation on a degraded image, and then obtaining a normalized spatial spectrum of the degraded image; reconstructing a frequency spectrum of an original image of the degraded image according to a spectral distribution rule of a natural image and estimating a normalized spatial spectrum of the original image; comparing the normalized spatial spectrum of the degraded image and a reconstructed frequency spectrum, and obtaining an optical transfer function of a system so as to obtain a point spread function of the degraded image; and carrying out Wiener filtering restoration on the degraded image according to the obtained point spread function so as to obtain an ideal image. The real-time blind image restoration method does not contain time-consuming iteration, and is simple and easy to implement, so that algorithm complexity is reduced, meanwhile, accurate estimation of the point spread function is obtained, and a good image restoration effect is realized.

Description

A kind of realtime graphic blind restoration method

Technical field

The invention belongs to technical field of image processing, relate to adaptive optical image and restore theory and means, be specifically related to the algorithm of blindly restoring image real-time in a kind of adaptive optics.

Background technology

Image restoration is the reason based on image degradation, utilizes the degraded image observed to recover real scene clearly.Image restoration technology is widely used in astronomical sight, remote sensing, many fields such as security monitoring and medical imaging.In most cases, the point spread function PSF (Point Spread Function) of image degradation is unknown, therefore carries out restoration disposal for degraded image and usually adopts blind restoration algorithm.Current blind restoration algorithm is iteration in essence mostly, due to the problem of some pathosis of degraded image data, often there will be local convergence or Divergent Phenomenon in the blind recuperation of iteration; Even if iterative process is stablized, iteration is also up to thousands of times even thousands of times, very consuming time, is difficult to the requirement meeting realtime graphic recovery.

When without any priori, it is very difficult for estimating PSF, until today, how to estimate that PSF is still the topic being active in this field exactly.Therefore, propose a kind of blind restoration method of realtime graphic fast and effectively and become current this area important technological problems urgently to be resolved hurrily.And the image degradation problem that G class point spread function causes is a kind of reason of common image degradation, because G class point spread function image degradation system occupy very large ratio, such as some astronomical sights, satellite sounding, aerospace detection, Electron micrographs and magnetic resonance image (MRI) system etc.

APEX algorithm goes out the PSF of system according to the feature assessment of the spectrum information of degraded image, then adopts SECB (Slow Evolution of Continuation Boundary, the slow evolution of continuum boundary) restored method to restore degraded image.APEX algorithm principle is simple, does not have iteration link, can restore some degraded images the short time, realtime graphic restoration disposal has the unexistent advantage of the blind restoration algorithm of conventional iterative.But APEX algorithm still exists some problems: first, an APEX algorithm default constant value replaces the spectrum signature of target image, and in the physical sense and unreasonable, therefore evaluated error is larger; Secondly, some parameters of APEX algorithm need artificially to determine, and these parameters must rule of thumb be selected, and usually need trial repeatedly just can obtain suitable value; Again, APEX method uses SECB method to rebuild image, and the method is easily subject to the impact of noise.APEX algorithm is when estimating the frequency spectrum of true picture, and adopt simple constant to substitute, this does not meet the spectral nature of true picture, causes larger evaluated error, affects recovery effects.Although people have carried out some improvement afterwards, do not tackle the problem at its root.In addition, when for astronomical degraded image and micro-degraded image (because the target and background contrast of these images is very large, and background is simple), APEX algorithm can carry out image restoration effectively; But when processing ground image, due to the complexity of background, restoration result is usually unsatisfactory.When estimating the frequency of true picture, these methods are all adopt simple constant to substitute, and this does not meet the spectral nature of true picture.

Summary of the invention

The technical problem to be solved in the present invention is: for overcoming iteration blind restoration method poor real in image restoration technology application, non-iterative blind restoration method estimates inaccurate shortcoming to point spread function, the present invention proposes the method for blindly restoring image based on image spectrum and frequency double-log relation.This method is improved based on the ultimate principle of APEX method, the normalized spatial spectrum of the original image of degraded image is estimated, tried to achieve point spread function is averaged on four straight lines, improve the accuracy of estimation, it is fast that this method restores speed, point spread function is estimated accurately, to achieve good recovery effect.

To achieve these goals, the present invention is realized by following technical proposals:

A kind of realtime graphic blind restoration method, comprises the following steps:

1) Fourier transform is carried out to degraded image, then obtain the normalized spatial spectrum of degraded image;

2) according to the spectrum distribution rule of natural image, the frequency spectrum of the original image of degraded image is rebuild, estimate the normalized spatial spectrum of original image;

3) normalized spatial spectrum of degraded image and reconstructed spectrum are compared, ask for the optical transfer function of system, and then ask for point spread function;

4) according to obtained point spread function, Wiener filtering recovery is carried out to degraded image.

Further, described step 1) in, Fourier transform is carried out to degraded image g (x) and obtains G (μ, v), G (μ, v) is normalized: G ^*(μ, v)=G (μ, v)/σ (1)

Wherein, G ^*the normalized frequency that (μ, v) is degraded image, G (μ, v) is blurred picture, and σ is normaliztion constant, gets σ=G 0,0, and namely the Fourier transform of degraded image is in the value of initial point.

Further, described step 2) in, for μ=0, the reconstruction formula of original image normalized frequency is:

$\ln \left|F^{*}(0,v)\right|=\left\{\begin{array}{cc}\ln \left|G^{*}(0,v)\right|& \left|v\right|<n\\ g_1+k^{*}\log \left|n\right|-k^{*}\log |v& n\&le;\sqrt{v}\&le;N-1\\ g_2+k^{*}\log \left|n\right|-k^{*}\log \left|v\right|& -N\&le;\sqrt{v}\&le;-n\end{array}\right.---\left(2\right)$

Wherein, k is the slope of the picture rich in detail estimated, G ^*the normalized spatial spectrum that (μ, v) is degraded image, F ^*(μ, v) original image normalized spatial spectrum for rebuilding, g ₁and g ₂for ln|G ^*(μ, v) | the value at place, * is convolution symbol, μ and v is the coordinate of domain space, and N is image size.

G (x) is carried out Fourier transform and obtains G (μ, v), G (μ, v) is normalized, obtains G ^*the normalized frequency that (μ, v) is degraded image.

Further, described g ₁, g ₂get ln|G ^*(0, n) |, ln|G ^*(0 ,-n) | at ln|G ^*(0, v) | on weighting smooth value.

Further, described step 3) in, the point spread function of degraded image is estimated, by following process implementation:

1) when the point spread function of system is G class point spread function, the frequency domain representation of G class point spread function is:

H(ξ，η)＝exp[α(μ ²+v ²) ^β]，α＞0，0＜β≤1， (3)

In formula, α and β is the parameter in G class point spread function, μ and v is the coordinate of domain space;

Image degradation system can be expressed as

G(μ，v)＝F(μ，v)H(μ，v)+N(μ，v)＝exp[-α(μ ²+v ²) ^β]F(μ，v)+N(μ，v) (4)

In formula, the optical transfer function that H (μ, v) is system, the noise that N (μ, v) is system;

Take the logarithm to after the degeneration system normalization of image:

ln|G ^*(μ，v)|＝ln|exp[-α(μ ²+v ²) ^β|F ^*(μ，v)+N ^*(μ，v)| (5)

Wherein, the noise N after normalization is ignored ^*(μ, v), then obtain following formula:

ln|G ^*(μ，v)|＝-α(μ ²+v ²) ^β]+ln|F ^*(μ，v)| (6)

(2) formula is substituted in (6), obtains-α | v| ^{2 β}value;

Wherein, α, β are the parameter value of the point spread function that will estimate;

2) least square method is adopted to carry out matching-α | v| ^{2 β}, the value of α and β can be obtained from the result of matching, thus obtain the point spread function causing image degradation.

Further, pass through to select μ=0 to the value of α and β, v=0, μ=v, μ=-v four rectilinear direction is estimated respectively, arithmetic average is asked for α, what embody due to β is the index characteristic of point spread function, so ask for geometrical mean to β, reasonable α and β value can be obtained like this, thus avoid the possibility of distortion estimator in a single direction.

Further, described step 4) in, according to obtained point spread function, carry out Wiener filtering recovery to degraded image, adopt least mean-square error filter method to carry out image reconstruction, its computing formula is:

$F(\mathrm{\&mu;},v)=\left[\frac{1}{H(\mathrm{\&mu;},v)}\frac{{\left|H(\mathrm{\&mu;},v)\right|}^2}{{\left|H(\mathrm{\&mu;},v)\right|}^2+S}\right]G(\mathrm{\&mu;},v)---\left(7\right)$

Wherein, S is the signal to noise ratio (S/N ratio) of blurred picture, F (μ, v) be the Fourier transform of original image, the Fourier transform that G (μ, v) is degraded image, H (μ, v) be the point spread function of system, namely the frequency of point spread function is expressed, μ and v is the coordinate of domain space.

This method is according to the logarithmic relationship of natural image, rebuild the original image frequency spectrum of degraded image, then for the relation of rebuilding image and degraded image, on image spectrum, point spread function is asked in four straight line weightings, restores degraded image according to tried to achieve point spread function.This method does not comprise iteration consuming time, and method is simple, is easy to realize, thus reduces algorithm complex, obtains point spread function accurately simultaneously and estimates, achieve good image restoration effect.

Its feature of the present invention is:

1) this algorithm utilizes the frequency spectrum of true picture and the double-log linear relationship of frequency to rebuild the frequency spectrum of original image, more has physical significance, more accurately than traditional algorithm.

2) this method is improved based on the ultimate principle of APEX method, estimates the normalized spatial spectrum of the original image of degraded image, averages, improve the stability of estimation to tried to achieve point spread function on four straight lines.

3) this method adopts Wiener filtering to replace SECB algorithm to restore degraded image, decrease the interactive session in recuperation, improve recovery speed, simultaneously, Wiener Filter Method is better to the inhibition of noise than SECB algorithm, improves the recovery effect of image.

Accompanying drawing explanation

Fig. 1 is the implementation method FB(flow block) that Real-time blind image of the present invention restores.

Fig. 2 rebuilds picture rich in detail spectrogram by with degraded image frequency spectrum.

Fig. 3 is the-α degeneration system that the original image frequency spectrum of reconstruction substitutes into image obtained | v| ^{2 β}the image of value.

Fig. 4 is p-α | v| ^{2 β}carry out the image that matching obtains α and β value.

Embodiment

Below in conjunction with drawings and Examples, further describing is done to the present invention.

As shown in Figure 1, realtime graphic blind restoration method of the present invention, comprises the following steps:

1) degraded image g (x) is carried out Fourier transform and obtains blurred picture G (μ, v),

Get σ=G (0,0), G (μ, v) is normalized,

G ^*(μ，v)＝G(μ，v)/σ (1)

Wherein, G ^*the normalized frequency that (μ, v) is degraded image, G (μ, v) is blurred picture, and σ is normaliztion constant, gets σ=F (0,0) ≈ G (0,0), and namely the Fourier transform of degraded image is in the value of initial point.

2) according to the spectrum distribution rule of natural image, rebuild the frequency spectrum of the original image of degraded image, the reconstruction formula (for μ=0) of original image frequency spectrum is

$\ln \left|F^{*}(0,v)\right|=\left\{\begin{array}{cc}\ln \left|G^{*}(0,v)\right|& \left|v\right|<n\\ g_1+k^{*}\log \left|n\right|-k^{*}\log |v& n\&le;\sqrt{v}\&le;N-1\\ g_2+k^{*}\log \left|n\right|-k^{*}\log \left|v\right|& -N\&le;\sqrt{v}\&le;-n\end{array}\right.---\left(2\right)$

Wherein, k is the slope of the picture rich in detail estimated, and k value is 0.8 to 1.2, n values when being 6 to 10, and the estimation effect obtained is better.G ^*the normalized spatial spectrum that (μ, v) is degraded image, F ^*(μ, v) original image normalized spatial spectrum for rebuilding.G ₁and g ₂for ln|G ^*(μ, v) | the value at place, * is convolution symbol, μ and v is the coordinate of domain space, and the value of n is 6-10, N is image size.

Meanwhile, due to ln|G ^*(0, v) | be concussion, so g ₁, g ₂get ln|G ^*(0, n) |, ln|G ^*(0 ,-n) | at ln|G ^*(0, v) | on weighting smooth value, the picture rich in detail frequency spectrum ln|F of reconstruction ^*(0, v) | as shown in Figure 2, in figure, asterisk point in top is the frequency spectrum of the clear picture rebuild, and bottom discrete point is the frequency spectrum of degraded image.

3) normalized spatial spectrum of degraded image and reconstructed spectrum are compared, ask for the optical transfer function of system, and then ask for the point spread function of system.

For μ=0, details are as follows for the process of asking for of point spread function:

1., when the point spread function of system is G class point spread function, the frequency domain representation of G class point spread function is:

H(ξ，η)＝exp[α(μ ²+v ²) ^β] (3)

In formula, α and β is the parameter in G class point spread function, μ and v is the coordinate of domain space;

By selecting the value of suitable α and β, H (ξ, η) can meet a large class optical transfer function, wherein corresponding gaussian density during β=1, and in seabed imaging, there is important application the aspects such as computer tomography; β=5/6 correspondence long exposure atmospheric turbulence is fuzzy, has important application in the many-side such as astronomical sight, remote sensing; The corresponding Lorentzian density function in β=1/2, has been applied in medical imaging and has carried out modeling for X ray scattering phenomenon.In a word, the β meeting 0 < β≤1 meets the feature of many electronics equipments.

Then the degeneration system of image is expressed as:

G(μ，v)＝F(μ，v)H(μ，v)+N(μ，v)＝exp[-α(μ ²+v ²) ^β]F(μ，v)+N(μ，v) (4)

In formula, the optical transfer function that H (μ, v) is system, the noise that N (μ, v) is system;

Take the logarithm to after the degeneration system normalization of image:

ln|G ^*(μ，v)|＝ln|exp[-α(μ ²+v ²) ^β]F ^*(μ，v)+N ^*(μ，v)| (5)

Wherein, the noise N after normalization is ignored ^*(μ, v), then obtain following formula:

ln|G ^*(μ，v)|＝-α(μ ²+v ²) ^β]+ln|F ^*(μ，v)| (6)

By ln|F ^*(0, v) | substitute into ln|G ^*(0, v) | ≈-α (v) ^{2 β}-ln|F ^*(0, v) |, can-α be obtained | v| ^{2 β}value, result of calculation is as shown in Fig. 2.

Wherein, α, β are the parameter value of the point spread function that will estimate;

2. least square method is adopted to carry out matching α | v| ^{2 β}, the value of α and β can be obtained from the result of matching, thus obtain the point spread function causing image degradation.

Due to-α after β is fixing | v| ^{2 β}bell, in-N≤v≤0 curve monotone increasing, in 0≤v≤N inner curve monotone decreasing.Therefore, the not all data of Fig. 3 can be used for the parameter of estimation point spread function, the most effective when only having the data estimation point spread function near initial point.Set a threshold value scale herein and scale < N, by the data in v ∈ [scale ,-scale] scope from-α | v| ^{2 β}in be extracted (as shown in Figure 4), scale is required that it makes got point all within the scope of two " troughs ".Discrete point is the-α that scale intercepts | v| ^{2 β}value, curve is-α | v| ^{2 β}matching.

After obtaining suitable data, adopt least square method to carry out matching-α | v| ^{2 β}, the value of α and β can be obtained from the result of matching, thus obtain causing image degradation point spread function.

In the degraded image of reality, due to the impact of degeneration factor and other enchancement factors, the spectrum information of the Turbulence-degraded Images of a frame reality can not Ω=(μ, v) | μ ²+ v ²≤ ω ²meet radiation symmetric condition completely in scope, can only approach to a certain extent or approximately meet radiation symmetric.Therefore, Ω plane was chosen the different straight lines of initial point, it is different for carrying out to log spectrum curve the value that matching draws, thus result in the instability that point spread function estimates.When the direction of selected straight line is larger by enchancement factor (as system noise, truncation error etc.) impact, likely causes the inefficacy of APEX method, be difficult to estimate rational α and β value.By selecting μ=0, v=0, μ=v, μ=-v four rectilinear directions are estimated respectively, ask for arithmetic average to α, what embody due to β is the index characteristic of point spread function, so ask for geometrical mean to β, reasonable α and β value can be obtained like this, thus avoid the possibility of distortion estimator in a single direction.

4) image restoration, the present invention adopts Wiener filtering to restore.

Apex algorithm adopts SECB method to carry out deconvolution image restoration usually.SECB method proposes for the feature of G class point spread function limitlessly detachable, but carries out interactive adjustment to parameter due to needs, relatively consuming time.The present invention adopts " least mean-square error filtering (Wiener filtering) " method to carry out image reconstruction.Wiener filtering obtains Recovery image by the square error minimized between original image and restored image, and can regard a special case of SECB method as, have stronger denoising effect, its computing formula is:

$F(\mathrm{\&mu;},v)=\left[\frac{1}{H(\mathrm{\&mu;},v)}\frac{{\left|H(\mathrm{\&mu;},v)\right|}^2}{{\left|H(\mathrm{\&mu;},v)\right|}^2+S}\right]G(\mathrm{\&mu;},v)---\left(7\right)$

Wherein, S is the signal to noise ratio (S/N ratio) of blurred picture, F (μ, v) be the Fourier transform of original image, the Fourier transform that G (μ, v) is degraded image, the optical transfer function that H (μ, v) is system, μ and v is the coordinate of domain space.

Claims (7)

1. a realtime graphic blind restoration method, is characterized in that, comprises the following steps:

1) Fourier transform is carried out to degraded image, then obtain the normalized spatial spectrum of degraded image;

2) according to the spectrum distribution rule of natural image, the frequency spectrum of the original image of degraded image is rebuild, estimate the normalized spatial spectrum of original image;

3) normalized spatial spectrum of degraded image and the ideal image normalized spatial spectrum of reconstruction are compared, ask for the optical transfer function of system, and then ask for the point spread function of degraded image;

4) according to obtained point spread function, Wiener filtering recovery is carried out to degraded image.

2. realtime graphic blind restoration method according to claim 1, is characterized in that, step 1) in, Fourier transform is carried out to degraded image g (x) and obtains G (μ, v), G (μ, v) is normalized:

G ^*(μ,v)＝G(μ,v)/σ (1)

Wherein, G ^*the normalized frequency that (μ, v) is degraded image, G (μ, v) is blurred picture, and σ is normaliztion constant, gets σ=F (0,0) ≈ G (0,0), and namely the Fourier transform of degraded image is in the value of initial point.

3. realtime graphic blind restoration method according to claim 1, is characterized in that, step 2) in, the reconstruction formula of original image normalized spatial spectrum is:

$\ln \left|F^{*}(\mathrm{\&mu;},v)\right|=\left\{\begin{array}{c}\ln \left|G^{*}(\mathrm{\&mu;},v)\right|\left|\sqrt{{\mathrm{\&mu;}}^2+v^2}\right|<n\\ g_1+k*\log \left|n\right|-k*\log \left|\sqrt{{\mathrm{\&mu;}}^2+v^2}\right|n\&le;\sqrt{{\mathrm{\&mu;}}^2+v^2}\&le;N-1\\ g_2+k*\log \left|n\right|-k*\log \left|\sqrt{{\mathrm{\&mu;}}^2+v^2}\right|-N\&le;\sqrt{{\mathrm{\&mu;}}^2+v^2}\&le;-n\end{array}\right.---\left(2\right)$

Wherein, k is the slope of the picture rich in detail estimated, G ^*the normalized spatial spectrum that (μ, ν) is degraded image, F ^*(μ, v) original image normalized spatial spectrum for rebuilding, g ₁and g ₂for ln|G ^*(μ, ν) | the value at place, * is convolution symbol, μ and ν is the coordinate of domain space, and N is image size.

4. realtime graphic blind restoration method according to claim 3, is characterized in that, described g ₁, g ₂get ln|G ^*(0, n) |, ln|G ^*(0 ,-n) | at ln|G ^*(0, v) | on weighting smooth value.

5. the realtime graphic blind restoration method according to any one of 1-3 claim, is characterized in that, step 3) in, the point spread function of degraded image is by following process implementation:

1) when the point spread function of system is G class point spread function, the frequency domain representation of G class point spread function is:

H(ξ,η)＝exp[α(μ ²+v ²) ^β],α＞0,0＜β≤1， (3)

In formula, α and β is the parameter in G class point spread function, μ and v is the coordinate of domain space;

Then the degeneration system of image is expressed as:

G(μ,ν)＝F(μ,ν)H(μ,ν)+N(μ,ν)＝exp[-α(μ ²+ν ²) ^β]F(μ,ν)+N(μ,ν) (4)

In formula, the optical transfer function that H (μ, v) is system, the noise that N (μ, ν) is system;

Take the logarithm to after the degeneration system normalization of image:

ln|G ^*(μ,v)|＝ln|exp[-α(μ ²+v ²) ^β]F ^*(μ,v)+N ^*(μ,v)| (5)

Wherein, the noise N after normalization is ignored ^*(μ, v), then obtain following formula:

ln|G ^*(μ,v)|＝-α(μ ²+v ²) ^β]+ln|F ^*(μ,v)| (6)

(2) formula is substituted in (6), obtains-α | v| ^{2 β}value;

Wherein, α, β are the parameter value of the point spread function that will estimate;

2) least square method is adopted to carry out matching-α | v| ^{2 β}, the value of α and β can be obtained from the result of matching, thus obtain the point spread function causing image degradation.

6. realtime graphic blind restoration method according to claim 5, is characterized in that, passes through to select μ=0, v=0, μ=v, μ=-v four rectilinear directions are estimated respectively, asks for arithmetic average, ask for geometrical mean to β to α the value of α and β.

7. realtime graphic blind restoration method according to claim 1, is characterized in that, step 4) in, adopt least mean-square error filter method to carry out image reconstruction, its computing formula is:

$F(\mathrm{\&mu;},v)=\left[\frac{1}{H(\mathrm{\&mu;},v)}\frac{{\left|H(\mathrm{\&mu;},v)\right|}^2}{{\left|H(\mathrm{\&mu;},v)\right|}^2+S}\right]G(\mathrm{\&mu;},v)---\left(7\right)$

Wherein, S is the signal to noise ratio (S/N ratio) of blurred picture, F (μ, v) be the Fourier transform of original image, the Fourier transform that G (μ, v) is degraded image, H (μ, v) be the optical transfer function of system, namely the frequency of point spread function is expressed, μ and v is the coordinate of domain space.