CN101877121A - Intermediate frequency based method for blindly restoring image - Google Patents

Abstract

The invention discloses a intermediate frequency based method for blindly restoring an image, comprising a concept of an intermediate frequency domain of the image and a method for geometrically positioning according to the characteristics. The degraded PSF (Payload-Structure-Fuel) can be quickly and accurately estimated by using the intermediate frequency domain and further the image can be quickly restored by using Wiener filter based on the intermediate frequency domain. Aiming at G type PSF suitable for a large amount of imaging systems, an image blind restoration method based on the intermediate frequency for estimating the G type PSF can be realized by using the technology and a single parameter fitting method. The reliability and the speed can be improved by using a special smooth and quick algorithm. The method has the characteristics of top speed, stability, fewer parameters and wide application, can carrying out blind processings, such as degradation resistance, noise suppression, detail enhancement and the like, on the image in real time and can be widely applied to the fields of astronomy, military affairs, medicines, remote sensing, televisions and the like.

Description

Method for blindly restoring image based on intermediate frequency

Technical field

The invention belongs to the blindly restoring image technical field, be specifically related to a plurality of technology neighborhoods such as the theory of realtime graphic processing, Deconvolution Technique, image detail enhancing, blurred picture sharpening, squelch etc. and algorithm.

Background technology

Image restoration is a hot spot technology in recent years, has been widely used for fields such as astronomy, military affairs, medical science, remote sensing, TV.Its purpose is to recover original clear scene from the observed image that degrades.The model that degrades on the mathematics can be described as original image and point spread function (point spread function, PSF) in the convolution that has under the situation of additive noise, as shown in the formula:

$g(x,y)=f(x,y)\⊗h(x,y)+n(x,y)---\left(1\right)$

X wherein, y is the two-dimensional coordinate on plane, g (x y) is the gray scale of degraded image, f (x, the y) gray scale of original picture rich in detail, h (x y) is the gray scale of PSF of degrading, n (x y) is additive noise, Represent convolution algorithm.By Fourier transform, this model can be at frequency domain representation

G(u，v)＝F(u，v)·H(u，v)+N(u，v)(2)

U wherein, v represents discrete frequency, G (u, v), F (u, v), H (u, v) and N (u, v) be respectively g (x, y), f (x, y), h (x, y) and n (x, Fourier transform y).Just (x, y) middle (x y) recovers f from g in image restoration.Yet, in actual conditions, h (x, y) and n (x, y) unknown often, so that (this situation promptly is called as blindly restoring image for x, y) unusual difficulty to recover f.

So far, people have proposed many blindly restoring image technology, such as zero page face partition method, the blind deconvolution method of iteration, simulated annealing, constraint recurrence liftering method, ARMA parameter estimation method, total least square method, the overall variational method, anisotropic diffusion equation method, minimum entropy deconvolution method or the like.Yet, use iteration or recursion step in these methods mostly, cause the possibility of result not restrained, speed is very slow, even may make picture quality worse and worse.Therefore, they almost can't be applied to needs to handle in real time in the engineering of image.

In order to realize the blind recovery of realtime graphic, people such as Carasso have proposed that a kind of (fastFourier transform, method FFT) are called the APEX method based on fast fourier transform.The APEX method can be finished the recovery of a width of cloth 512*512 image in tens of seconds, and good effect is arranged.Also occur the improvement and the application of some these methods in recent years, also obtained certain success.The APEX method improves the speed major reason and is that it is set at a kind of function of very extensively seeing with PSF, is called the G class function, and its Fourier transform has following form:

$H(u,v)={\&Integral;}_{R^2}h(x,y)\exp [-\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="H2009102359465E0000071.png,H2009102359465E0000072.png"\; state="\{\{state\}\}">i</figure-callout>}2\mathrm{\&pi;}(\mathrm{xu}+\mathrm{vy})]\mathrm{dxdy}---\left(3\right)$

$=\exp [-c{({\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="H2009102359465E0000071.png,H2009102359465E0000072.png"\; state="\{\{state\}\}">u</figure-callout>}}^2+v^2)}^{\mathrm{\&beta;}}],c>\mathrm{0,0}<\mathrm{\&beta;}\&le;1$

Wherein, β is relevant with imaging circumstances, has determined fuzzy kind, and c is relevant with imaging process, has determined fuzzy degree.Very meaningfully, this function appears in a large amount of physical processes, it has not only comprised common Gaussian process (β=1) and Lorentzian process (β=1/2), can also portray the optical transfer function of many imaging systems, as: the X ray in the also corresponding actinoscopy in β=1/2 is scattering into picture, the long exposure image of the corresponding atmospheric turbulence in β=5/6, the optical transfer function of the approximate corresponding ideal lens limited diffraction in β=3/4, multiple optoelectronic device can be portrayed in 1/2≤β≤1, and nearly all photographic film all has the characteristic of 0＜β≤1, or the like.The APEX method is the parameter c by estimating G class PSF and the β purpose that reaches blindly restoring image just.

Yet still there are many problems in the APEX method.At first, the APEX method has too much artificial input parameter, and these important parameters can only be selected by experience, need repeatedly interactively trial just can obtain suitable value (and not necessarily can find) sometimes, even may make problem become difficult further.Secondly, the APEX method when estimated parameter with the ln|F of the unknown (u, v) | directly replace, and (least squaresestimation LSE) asks c and β with nonlinear least square method with a constant.But because ln|F (u, v) | and ln|G (u, v) | generally all is the concave curve that the overall situation singly subtracts, thus the β that obtains always less than 0.5, this often and actual conditions inconsistent.Though also may realize to a certain degree recovery with c that obtains and β, this has greatly limited its range of application.Once more, among the non-linear LSE recursive procedure is arranged, this will reduce the speed of restoring.At last, the APEX method is used the slow evolution of continuum boundary (slow evolution of continuation boundary, SECB) technology is carried out the deconvolution reconstructed image, this technology must be attempted a series of time value and be selected optimum, and each trial all needed about one second, had wasted many times.Though people improved to some extent afterwards, do not tackle the problem at its root, more consuming time on the contrary.Therefore, in the stronger application of real-time, these methods still can not meet the demands.

Blindly restoring image except h (x, y) and n (x y) outside the Wei Zhi difficult point, also has a difficulty to be that it is the inverse problem of a deconvolution, i.e. the slight disturbance of noise all may cause restored image to distort greatly.So also having proposed the technology of many deconvolution, people overcome this point, such as SECB technology above-mentioned.But the SECB technology proposes at the characteristics of G class function limitlessly detachable specially, and therefore certain limitation is arranged.Having a kind of Deconvolution Technique of classics to be called in addition " least mean-square error filtering (Wiener filtering) ", is to derive by the square error that minimizes between original image and the restored image, and its restored image calculates with following formula

$\hat{F}(u,v)=\frac{{\hat{H}}^{*}(u,v)G(u,v)}{{\left|\hat{H}(u,v)\right|}^2+\&prod;(u,v)}---\left(4\right)$

Wherein

With

Be respectively the Fourier transform estimation of restored image and PSF,

Be

Complex conjugate, ∏ (u, v)=| N (u, v) | ²/ | F (u, v) | ²Be the ratio of noise and ideal image power spectrum.(u when v) unknown, can replace with a constant Γ, and promptly following formula is approximately as ∏

$\hat{F}(u,v)\&ap;\frac{{\hat{H}}^{*}(u,v)G(u,v)}{{\left|\hat{H}(u,v)\right|}^2+\mathrm{\&Gamma;}}---\left(5\right)$

Wherein Γ is the ratio of noise to the spectral density of signal, in the time of statistical property the unknown of noise, and the key of determining just to become decision recovery effect quality of Γ value.People have proposed the optimum value that certain methods is determined Γ to this, also obtain certain achievement.It is exactly that operand is bigger that but existing method all has a common shortcoming, need calculate Γ value with the statistical information or the recursive operation of image, so has reduced the efficient of Wiener filtering on the contrary, and this also is this technology part that haves much room for improvement.Yet, though whether optimum does not still have final conclusion in Wiener filtering, it simply fast advantage make it to come compared with other technologies, still be more suitable in real-time application.

Summary of the invention

Technology of the present invention is dealt with problems: overcome the deficiency that the blind recovery technique speed of conventional images is slow excessively, range of application is narrow, a kind of method for blindly restoring image based on intermediate frequency is provided, this method makes the real-time application of blindly restoring image become possibility, also has the removal of images noise and strengthen in the image function such as point target in quick restored image.

Technical scheme of the present invention: based on the method for blindly restoring image of intermediate frequency, performing step is as follows:

(1) frequency domain that contains PSF information in the image spectrum is positioned, described frequency domain localization method is:

(1.1) FFT and the standardization of calculating blurred picture;

(1.2) some of calculating FFT are crossed the frequency domain separation of former dotted line, thereby determine the scope of frequency domain;

(2) utilize frequency domain to estimate G class PSF, method is:

(2.1) according to the value of model index β in the environmental selection G class function;

(2.2) frequency domain of the identical former dotted line of mistake is carried out smoothly, and estimate the derivative of PSF in the frequency domain with the linear LSE of one-parameter;

(2.3) estimated value of fog-level parameter c in the usefulness derivative calculations G class function of PSF

Obtain the estimator of G class PSF;

(3) utilize frequency domain to optimize in the Wiener filtering noise spectral density than the value of Γ, and with Wiener filtering restored image.

Described step (1.2) calculate FFT some cross being achieved as follows of frequency domain separation of former dotted line:

(1.2.1) logarithm of crossing former dotted line is carried out " moving right average " smoothly, promptly the value of each point is replaced with the mean value of its several consecutive point of the right, with the straight line connection end point and calculate the point of this straight line ultimate range, promptly obtain the separation between frequency domain and the high-frequency domain then;

(1.2.2) " mobile left mean " carried out smoothly in the medium and low frequency territory of crossing former dotted line, promptly the value of each point is replaced with the mean value of its several consecutive point of left side, with the straight line connection end point and calculate the point of this straight line ultimate range, promptly obtain the separation between lower frequency region and the frequency domain then;

(1.2.3) zone between two separations calculating of above step is exactly this frequency domain of crossing former dotted line, when the PSF invariable rotary, then with after the frequency domain annulusization of this line as the frequency domain of image, otherwise then determine (when getting two former dotted lines of mistake, frequency domain is an elliptical ring) jointly by many former dotted lines of mistake.

Utilize frequency domain to optimize noise spectral density in the Wiener filtering in the described step (3) to be than the value method of Γ:

(3.1) utilize the middle and high frequency domain separation of the identical former dotted line of mistake to make the adaptive best Wiener filtering of Γ;

(3.2) introduce squelch parameter η, with finely tune Γ to be in harmonious proportion the contradiction of noise and image detail: η is comparatively clear near 1 o'clock noise and details, noise and details weakened simultaneously to some extent when η reduced, and therefrom chose an equilibrium value and restored with the best that realizes image.

The present invention's beneficial effect compared with prior art is:

(1) on the speed, since without any recurrence and iterative process, also without any matrix operation, even has only few power operation, therefore the time of restoring almost only consumes in necessary FFT and Wiener filtering, make recovery speed near the limit, than the blind recovery technique of conventional images more than fast ten times.Restore such as the image to a secondary 512*512, other prior art is the several seconds at least, several at most branches, and the present invention is at common computer (CPU:Intel Pentium 4,2.40GHz; Memory:512MB) 0.5 second time spent only on, and if with equipment such as DSP, recovery time even can reach Millisecond.

(2) on the accuracy, the proposition of " frequency domain " notion and the location of quantification thereof make the accurate estimation of PSF become possibility, here be the picture rich in detail of two width of cloth 512*512 with some examples explanation: Fig. 1 (a) and (d), at first using c=0.075 and β=0.5 that (a) carried out computer simulation degrades and adds to make an uproar and obtain (b), restore with β=0.5 and η=0.8 then and obtain (c), wherein calculate

Very near 0.075 of reality, this is the unapproachable accuracy of prior art; (d) carried out computer simulation with c=0.0025 and β=5/6 degrade that adding makes an uproar obtains (e), obtain (f) after restoring with β=5/6 and η=0.8 again, wherein calculate

Too very near actual value.Experiment shows that as long as model is suitable, the PSF that estimates will be very accurate, thereby obtains satisfied recovery effect.

Recovery speed and effect that it is good also are confirmed in the photo of actual photographed, as Fig. 2-5.First width of cloth is actual fuzzy photo among each figure, and second width of cloth is the picture after restoring.As seen, photographic quality obviously improves, and details is more clear, and acutance has improved nearly 100%; Recovery speed is also quite fast, and the image of a width of cloth 1280*1024 is consuming time not even above 2.5 seconds on common computer.Listed the detail statistics data of Fig. 2-5 in the table 1, wherein the processing time is added up on the same computer, and acutance uses the sharpness evaluation function based on the Roberts edge detection operator to calculate.

The statistics of table 1 Fig. 2-5

(3) on the simplicity because the present invention has only two input parameters, and they change among a small circle result's influence very little, so needn't as existing some method, a large amount of mutual test of needs come various parameter is chosen.Parameter beta among the present invention is used for adaptive to true environment, thereby has widened the range of application (nearly all blurred picture all have a suitable β to can be used to carry out sharpening handle) of algorithm; And even more advantageously, common same environment is fit to same β value, so as long as obtained suitable β value, just need not to adjust again.Another squelch parameter η is used to be in harmonious proportion the contradiction of noise and image detail, and η is comparatively clear near 1 o'clock noise and details, and noise is inhibited when η reduces, but details also has loss slightly.But it is also little to result's influence, is applicable to most applications about η=0.8 substantially, unless noise is too obvious, otherwise also need not to adjust.This also is one of advantage of this invention, can be general, and also fine-tuning.

(4) on the function because the adjustability of parameter beta, the present invention also have that other prior art is difficult to realize specific function, strengthen as squelch and point target.As long as the β that chooses is suitable, just can suppress ground unrest, thus outstanding wherein point target.Fig. 6 (a) is that a width of cloth is taken from the 1024*1024 of astro tracker picture, because The noise almost be can't see any target.But after using β=0.8 to handle, noise weakens significantly, and plural target becomes high-visible, as Fig. 6 (b).Fig. 6 (c) and (d) be respectively (a) and three-dimensional space grid figure (b), target is outstanding more obvious among Fig. 6 (d).

In a word, the characteristics that the present invention has is very fast, stable, parameter is few, purposes is wide can instead degrade to realtime graphic, blind processing such as squelch, details enhancing, can be widely used in every field such as astronomy, military affairs, medical science, remote sensing, TV.

The drawing explanation

Fig. 1 a, Fig. 1 b, Fig. 1 c, Fig. 1 d, Fig. 1 e, Fig. 1 f are the recovery effect figure of computer simulation degraded image;

Fig. 2 is the recovery effect figure of the fuzzy photograph of a width of cloth 440*330;

Fig. 3 is the recovery effect figure of the fuzzy photograph of a width of cloth 640*480;

Fig. 4 is the recovery effect figure of the fuzzy photograph of a width of cloth 1024*768;

Fig. 5 is the recovery effect figure of the fuzzy photograph of a width of cloth 1280*1024;

Fig. 6 a, Fig. 6 b, Fig. 6 c, Fig. 6 d are recovery effect figure and the space networks trrellis diagrams thereof that 1024*1024 contains the point target photograph;

Fig. 7 is the frequency domain curved surface comparison diagram of image spectrum and PSF frequency spectrum;

Fig. 8 is the localization method key diagram in vision intermediate frequency territory;

Fig. 9 is typical T ₀(u), T ₁(u) and T ₂(u) comparison diagram of curve;

Figure 10 the present invention is based on the process flow diagram that intermediate frequency is estimated the method for blindly restoring image of G class PSF.

Embodiment

The present invention comprises three steps, and its specific implementation is as follows:

1. the location of image " frequency domain ".

In general, the most energy of piece image all concentrates near the low frequency region its frequency spectrum initial point, and noise pollution the high-frequency region at its frequency spectrum edge.Therefore, directly utilize the FFT of image, promptly (u estimates that v) (u, (u v) almost is impossible to H thereby v) recover F to G.Yet (u, information v) and not polluted by noise are referred to as " frequency domain " of image here, as shown in Figure 7 to find to have between lower frequency region and high-frequency domain a special region to keep H.The zone of black is exactly the frequency domain R of image among Fig. 7 _M, it is rendered as a ring-type; The non-black regional centralized of the ring heart most of energy of image, i.e. the lower frequency region of image; The outer zone of ring is polluted by noise, i.e. the high-frequency domain of image.Image spectrum ln|G from Fig. 7 (u, v) | and PSF frequency spectrum ln|H (u, v) | contrast as seen, the two smooth surface in frequency domain RM is extremely similar, (u, frequency domain has v) kept H (u, raw information v) to G in other words.So,, just can realize the estimation of PSF more easily if can determine the scope of frequency domain quantitatively.The technical matters that first step of the present invention will solve is exactly how to determine the scope of basic, normal, high frequency domain, promptly frequency domain is positioned.The present invention proposes a new method of geometry that does not contain any artificial parameter, and the location that can realize frequency domain effectively is as follows:

If the standardization FFT conversion G of blurred picture (u v) is the complex matrix of a M * N, selects the mould value of its a former dotted line of mistake, such as | G (u, 0) |.Because G (u, 0) is the conjugation symmetry, so only need consider half period u ∈ D _T=[0, u _T] on value just passable, wherein

Expression rounds to zero.Lower frequency region is located near the u=0 so, and high-frequency domain is positioned at u=u _TNear.At | G (u, 0) | curve in, two tangible overall turning points are arranged usually, they are the frontier point of frequency domain just.On the one hand, because lower frequency region has been assembled most of energy, so | G (u, 0) | ²Or | G (u, 0) | curve in can clearly see an overall turning point, promptly low, frequency domain separation, its coordinate u _LMExpression; On the other hand,, and polluted by approximate equally distributed noise because the value of high-frequency domain is very little, so in log|G (u, 0) | curve in can find another overall turning point, promptly middle and high frequency domain separation, its coordinate u _MHExpression.So just can calculate u according to the turnover characteristic of curve _LMAnd u _MHWith the location frequency domain, Fig. 8 has presented basic ideas: at first use straight line (to use L _T(u) expression) connect ln|G (u, 0) | at D _TIn two end points, ln|G (u, 0) then | with L _T(u) form an inverted triangle, so ln|G (u, 0) | to L _T(u) point of ultimate range is exactly a vertex of a triangle, promptly overall turning point u _MHThen [0, u _MH] in to curve | G (u, 0) | carry out identical step, just can obtain u _LM

Yet, | G (u, 0) | normally vibration strongly of curve, considerable local turning point is arranged, so that can't correctly find out u _LMAnd u _MHOne of solution is a smooth curve, such as the method for moving average.But smoothly can redistribute energy, the result who leads to errors.In order to address this problem, the present invention uses a kind of special smoothing method, is called " moving right (left side) on average " method and calculates u _LMAnd u _MH, as follows:

At first consider u _MH, being penetrated into the high-frequency domain on its right for fear of the energy of frequency domain, right average out to is moved in definition:

$G_r\left(u\right)=\underset{\left(\mathrm{right\; average}\right)}{\mathrm{smooth}}\left[\ln \right|G(u,0)],u\&Element;D_T---\left(6\right)$

Promptly

$G_r\left(u\right)=\frac{1}{s_r}\underset{i=0}{\overset{s_r-1}{\mathrm{\&Sigma;}}}\ln \left|G(u+i,0)\right|,u\&Element;D_T---\left(7\right)$

S wherein _rBe to move right average span, the number of the consecutive point that promptly are averaged can be used D _TWidth determine, as

So just can use G _r(u) find out u _MH: at first cross (0, G _rAnd (u (0)) _T, G _r(u _T)) 2 make straight line L _T(u), G so _r(u) last is exactly u from this straight line point farthest _MHBut because G _r(u _T) polluted the most handy G by noise at high-frequency domain _r(u) average

Replace G _r(u _T).Therefore, L _T(u) satisfy

$[G_r\left(0\right)-{\stackrel{\&OverBar;}{G}}_r]u+u_T{L}_T\left(u\right)-u_T{G}_r\left(0\right)=0---\left(8\right)$

And G _r(u) point on is to L _T(u) distance is

$\left|\right[G_r\left(0\right)-{\stackrel{\&OverBar;}{G}}_r]u+u_T{G}_r\left(u\right)-u_T{G}_r\left(0\right)|/K_T---\left(9\right)$

Wherein

Be a constant, judge and to remove when being worth most.Like this, just can obtain

$u_{\mathrm{MH}}=\arg \underset{u\&Element;D_T}{\max }\left|\right[G_r\left(0\right)-{\stackrel{\&OverBar;}{G}}_r]u+u_T{G}_r\left(u\right)-u_T{G}_r\left(0\right)|---\left(10\right)$

Like this, even to G _r(u) carry out pointwise and calculate u _MHAlso can calculate very soon.

Then, because 0≤u _LM≤ u _MH, calculate u _LMThe time only need consider field of definition u ∈ D _LM=[0, u _MH].In order to make the energy on the lower frequency region of the left side be independent of frequency domain, define mobile left mean and be:

$G_l\left(u\right)=\underset{\left(\mathrm{left\; average}\right)}{\mathrm{smooth}}\left[\right|G(u,0)\left|\right],u\&Element;D_{\mathrm{LM}}---\left(11\right)$

Promptly

$G_l\left(u\right)=\frac{1}{s_l}\underset{i=0}{\overset{s_l-1}{\mathrm{\&Sigma;}}}\left|G(u-i,0)\right|,u\&Element;D_{\mathrm{LM}}---\left(12\right)$

S wherein ₁It is the span of mobile left mean.Equally, G _lDistance mistake (u) (0, G _lAnd (u (0)) _MH, G _r(u _MH)) straight line point farthest, be exactly low, frequency domain separation u _LM, promptly

$u_{\mathrm{LM}}=\arg \underset{u\&Element;D_{\mathrm{LM}}}{\max }\left|\right[G_l\left(0\right)-G_l\left(u_{\mathrm{MH}}\right)]u+u_{\mathrm{MH}}{G}_l\left(u\right)-u_{\mathrm{MH}}{G}_l\left(0\right)|---\left(13\right)$

Like this, the frequency domain that belongs to G (u, 0) is D just _M=[u _LM, u _MH].If | H (u, v) | be (such as the G class function) of invariable rotary, so | G (u, v) | also approximate invariable rotary, at this moment (u, frequency domain v) can be considered annulus to G Otherwise (u, the former dotted line of mistake v) calculates corresponding frequency domain, last common decision G (u, frequency domain v) need to choose many G.Such as get two former dotted line G of mistake (u, 0) and G (0, in the time of v), if calculate separately frequency domain for [u _LM, u _MH] and [v _LM, v _MH], (u, frequency domain v) then is an elliptical ring to G so

Deng.

Because frequency domain positions according to its specific modality, so said as the front, it not only can ignore The noise, and avoided F (u, thus most energy v) kept H (u, partial contour v) promptly have following characteristic:

|G(u，v)|＞＞|N(u，v)|，(u，v)∈R _M (14)

With

smooth[|G(u，v)|]≈smooth[C?|H(u，v)|]，(u，v)∈R _M (15)

Wherein C is a constant.Thereby, just can utilize these two character that PSF is estimated.

2. utilize frequency domain to estimate PSF (the present invention uses G class PSF).

Utilize frequency domain, can carry out blind estimation to PSF easily, and then realize the recovery of image.Because it is very extensive that G class PSF uses, and the characteristics of invariable rotary are arranged, be fit to very much estimate with the frequency domain method, so the present invention has proposed an one-parameter derivative fitting process fast to this, as follows:

At first (u v) chooses the former dotted line of identical mistake, is order for simplicity to H

H _a(u)＝ln|H(u，0)|＝-c(u ²) ^β(16)

Wherein c and β are for treating estimated parameter.Yet parameter beta is relevant with the environment that degrades, and (u, v) the profile of frequency domain is to be difficult to correct estimation only to utilize G.So the present invention is used as β an input parameter without hesitation, directly comes the people for choosing according to environment.The benefit of doing like this is, not only makes c become unique estimated parameter for the treatment of, thereby and can make β break away from restriction to mate truth better, widened range of application.More fortunately, nearly all blurred picture all has a suitable β to can be used to carry out the sharpening processing.According to (15) and (16) formula, after choosing β, can use

$G_a\left(u\right)=\underset{\left(\mathrm{average}\right)}{\mathrm{smooth}}\left[\ln \right|G(u,0)\left|\right],u\&Element;D_M---\left(17\right)$

Estimate corresponding c value.Here

Represent conventional moving average level and smooth, promptly

S wherein _aIt is the span (odd number) of moving average.In frequency domain, because noise can ignore, and ln|G (u, 0) |, ln|F (u, 0) | and H _a(u) after level and smooth, can be similar to, so have with straight line

k _Gu+b _G≈smooth[ln|G(u，0)|]

≈smooth[H _a(u)]+smooth[ln?|F(u，0)|](19)

≈k _Hu+b _H+k _Fu+b _F，u∈D _M

K wherein _G, k _HAnd k _FBe the slope of each straight line, b _G, b _HAnd b _FIt is the intercept of each straight line.The both sides differentiate gets

k _G≈k _H+k _F (20)

At first, because energy seldom and successively decrease in the frequency domain, so have

k _F＝sgn(u _LM-u _MH)·ε(21)

Wherein sgn () is a sign function, and ε is a very little positive number.Desirable for the sake of simplicity k _F=0, and its error can be eliminated when adjusting squelch parameter η.

Secondly, because H _a(u) singly subtract, so k _HMust equal H _a(u) at D _MIn the derivative of certain point.If use u _MRepresent this point, can determine its position like this: establishing has straight line in addition

Simultaneously by (u _LM, H _a(u _LM)) and (u _MH, H _a(u _MH)) 2 points, following relation is then arranged:

$H_a\left(u_{\mathrm{MH}}\right)=-c{\left(u_{\mathrm{MH}}^2\right)}^{\mathrm{\&beta;}}={\tilde{L}}_H\left(u_{\mathrm{MH}}\right)=k_H{u}_{\mathrm{MH}}+{\tilde{b}}_H---\left(22\right)$

$H_a\left(u_{\mathrm{LM}}\right)=-c{\left(u_{\mathrm{LM}}^2\right)}^{\mathrm{\&beta;}}={\tilde{L}}_H\left(u_{\mathrm{LM}}\right)=k_H{u}_{\mathrm{LM}}+{\tilde{b}}_H---\left(23\right)$

And

${\tilde{L}}_H^{\&prime;}\left(u_M\right)=k_H=H_a^{\&prime;}\left(u_M\right)=-2\mathrm{\&beta;c}\&CenterDot;\mathrm{sgn}\left(u_M\right)\&CenterDot;{\left|u_M\right|}^{2\mathrm{\&beta;}-1}---\left(24\right)$

Cancellation c, k _HWith

Can solve

$u_M=\mathrm{sgn}\left(u_M\right){[\frac{\mathrm{sgn}\left(u_M\right)}{2\mathrm{\&beta;}}\&CenterDot;\frac{{\left|u_{\mathrm{MH}}\right|}^{2\mathrm{\&beta;}}-{\left|u_{\mathrm{LM}}\right|}^{2\mathrm{\&beta;}}}{u_{\mathrm{MH}}-u_{\mathrm{LM}}}]}^{\frac{1}{2\mathrm{\&beta;}-1}},\mathrm{\&beta;}\&NotEqual;0.5---\left(25\right)$

When β=0.5, use u _M=(u _MH+ u _LM)/2 get final product.

In addition, slope k _GCan obtain with the linear LSE of one-parameter at an easy rate:

$k_G=\frac{\underset{i=u_{\mathrm{LM}}}{\overset{u_{\mathrm{MH}}}{\mathrm{\&Sigma;}}}[G_a\left(i\right)-G_a\left(u_M\right)](i-u_M)}{\underset{i=u_{\mathrm{LM}}}{\overset{u_{\mathrm{MH}}}{\mathrm{\&Sigma;}}}{(i-u_M)}^2}---\left(26\right)$

At last, with k _G, k _HAnd k _FThe estimated value that just obtains c behind the c is returned and solved to generation

For

$\hat{c}=\frac{-\mathrm{sgn}\left(u_M\right)}{2\mathrm{\&beta;}}\&CenterDot;{\left|u_M\right|}^{1-2\mathrm{\&beta;}}\&CenterDot;k_G---\left(27\right)$

Thereby produce H (u, being estimated as v)

$\hat{H}(u,v)=\exp [-\hat{c}{({\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="H2009102359465E0000071.png,H2009102359465E0000072.png"\; state="\{\{state\}\}">u</figure-callout>}}^2+v^2)}^{\mathrm{\&beta;}}]---\left(28\right)$

3. based on the Wiener filtering of frequency domain.

After estimating PSF, just can carry out image restoration with Wiener filtering.The parameter Γ of Wiener filtering is the key of recovery effect quality.By the location of frequency domain, the present invention proposes a method of simply and quickly determining the Γ value, and is as follows:

At first choose the former dotted line of identical mistake, and order

$T_0\left(u\right)=\ln [{\left|\hat{H}(u,0)\right|}^2+\mathrm{\&Pi;}(u,0)]---\left(29\right)$

$T_1\left(u\right)=\ln [{\left|\hat{H}(u,0)\right|}^2+\mathrm{\&Gamma;}]---\left(30\right)$

$T_2\left(u\right)=\ln \left[{\left|\hat{H}(u,0)\right|}^2\right]---\left(31\right)$

Fig. 9 is typical T ₀(u), T ₁(u) and T ₂(u) comparison diagram of curve, therefrom visible T ₀(u) and T ₁(u) only at high-frequency domain (with u=u _MHBe the boundary) bigger difference arranged.Like this, problem just becomes and how to choose Γ and make T ₁(u)) adaptive best T on the basis that suppresses noise ₀(u).Because D _LMMiddle noise can be ignored, so always have

T ₀(u)≈T ₁(u)≈T ₂(u)，u∈D _LM(32)

So u=u _MHIn time, have

T ₀(u _MH)≈T ₁(u _MH)≈T ₂(u _MH)(33)

Evenly distribute because noise is approximate, so at high-frequency domain D _H=(u _MH, u _T] in T is arranged ₀(u) ≈ T ₀(u _MH).Optimize Wiener filtering, then wish T ₁(u) ≈ T ₀(u) at D _HIn also can set up, but fact proved T ₁(u) 〉=T ₀(u) more help to suppress noise.What therefore, need is

T ₁(u)≥T ₀(u)≈T ₀(u _MH)≈T ₂(u _MH)，u∈D _H (34)

Following formula is set up, only needed

Or

$\underset{u\&Element;D_T}{\min }\left[T_1\left(u\right)\right]=T_2\left(u_{\mathrm{MH}}\right)+\mathrm{\&delta;}---\left(35\right)$

Get final product, wherein δ is a very little positive number.Simultaneously in order to guarantee

Order

$\mathrm{\&delta;}=(1-\mathrm{\&eta;})\left\{\underset{u\&Element;D_T}{\max }\right[T_2\left(u\right)]-T_2\left(u_{\mathrm{MH}}\right)\}---\left(36\right)$

Wherein η ∈ (0,1] be a squelch parameter.So obtain

$\underset{u\&Element;D_T}{\min }\left[T_1\left(u\right)\right]=T_2\left(u_{\mathrm{MH}}\right)+(1-\mathrm{\&eta;})\left\{\underset{u\&Element;D_T}{\max }\right[T_2\left(u\right)]-T_2\left(u_{\mathrm{MH}}\right)\}---\left(37\right)$

With T ₀(u) and T ₁(u) expression formula substitution also solves Γ and just obtains

$\mathrm{\&Gamma;}={\left|\hat{H}(u_{\mathrm{MH}},0)\right|}^{2\mathrm{\&eta;}}\&CenterDot;{\left[\underset{u\&Element;D_T}{\max }\right|\hat{H}(u,0)\left|\right]}^{2(1-\mathrm{\&eta;})}-{\left[\underset{u\&Element;D_T}{\min }\right|\hat{H}(u,0)\left|\right]}^2---\left(38\right)$

Like this, just utilize frequency domain to obtain the suitable Γ value of Wiener filtering, its calculating is very succinct, settles at one go.Squelch parameter η wherein is used to be in harmonious proportion the contradiction of noise and image detail, and η is comparatively clear near 1 o'clock noise and details, and noise is inhibited when η reduces, but details also has loss slightly.But it is little to result's also influence, is applicable to most applications about η=0.8 substantially, unless noise is too obvious, otherwise need not to adjust.When using G class PSF, because

So following formula can abbreviation be

$\mathrm{\&Gamma;}={\left|\hat{H}(u_{\mathrm{MH}},0)\right|}^{2\mathrm{\&eta;}}-{\left|\hat{H}(u_T,0)\right|}^2---\left(39\right)$

Again because

Be real function, so Wiener filtering is used

$\hat{F}(u,v)\&ap;\frac{G(u,v)}{\hat{H}(u,v)+\mathrm{\&Gamma;}/\hat{H}(u,v)}---\left(40\right)$

Form will be more quick.

By above three steps, blurred picture just can be restored fast.But notice to have three level and smooth computings in the above process, its speed will have a direct impact the performance of whole algorithm.If carry out smoothly with its definition or wave filter, the number of times of computing will be directly proportional with its span.Such as moving right average computation G _r(u) time, then to be s at least if directly calculate _r(u _T+ 1) sub-addition (is noted factor 1/s _rCan not take advantage of, because use G _r(u) calculate u _MHThe time can look like K _TRemove like that), work as s _rCalculated amount is with considerable when very big.Therefore, the present invention also proposes an additional quick and smooth method, and is as follows:

Calculate G _r(u) time, at first set up a new function I _r(u):

$\left\{\begin{array}{c}{I}_r(u_T+s_r)=0\\ I_r\left(u\right)=I_r(u+1)+\ln \left|G(u,0)\right|,u\&Element;[0,u_T+s_r-1]\end{array}\right.---\left(41\right)$

G then _r(u) just can calculate with following formula:

$G_r\left(u\right)=\frac{1}{s_r}[I_r\left(u\right)-I_r(u+s_r)],u\&Element;D_T---\left(42\right)$

Factor 1/s wherein _rCan save.Like this, to any s _rOnly need be (2u _T+ s _r+ 1) sub-addition (comprises and calculates I _r(u)) just can obtain G _r(u).While G _r(u) average

Also can directly calculate:

${\stackrel{\&OverBar;}{G}}_r=[I_r\left(0\right)-I_r(u_T+1)]/(u_T+1)---\left(43\right)$

In like manner, calculate G _l(u) also set up a new function I time earlier _l(u):

$\left\{\begin{array}{c}{I}_l(-s_l)=0\\ I_l\left(u\right)=I_l(u-1)+\left|G(u,0)\right|,u\&Element;[1-s_l,u_{\mathrm{MH}}]\end{array}\right.---\left(44\right)$

Like this to any s _lOnly need be (2u _MH+ s _l+ 1) sub-addition (comprises and calculates I _l(u)) just can obtain G _l(u) (can omit factor 1/s _l):

$G_l\left(u\right)=\frac{1}{s_l}[I_l\left(u\right)-I_l(u-s_l)],u\&Element;D_{\mathrm{LM}}---\left(45\right)$

To G _a(u), can be by I _r(u) calculate.To any s _a, only need be (u with following formula _MH-u _LM+ 1) sub-addition just can obtain G _a(u):

Factor 1/s wherein _aThough can not save, can be left to and calculate k _GDisposable being multiplied by in back.

Note span s _r, s _lAnd s _aAll the most handy level and smooth separately peak width decides.The present invention gets

With

Summary is got up, and the location of frequency domain and be a kind of basic theory based on the Wiener filtering of frequency domain is mainly used in the popularization of method; When it uses G class PSF, just form the blindly restoring image algorithm of estimating G class PSF based on intermediate frequency, can specifically be applied to every field.This algorithm can realize that also but hardware is realized by software, can realize devices realizations such as also available DSP, FPGA on computers.The flow process of this algorithm specifically can be carried out according to following steps as shown in figure 10:

1) according to occasion select β and η value (common same environment is fit to same β value, therefore can be earlier mode by trial obtain the best value, as if uncertain β=0.5 and η=0.8 of then advising);

2) FFT and the standardization of calculating degraded image;

3) use formula (47), (41), (42), (43) and (10) to calculate s respectively _r, I _r(u), G _r(u),

And u _MH

4) use formula (48), (44), (45) and (13) to calculate s respectively ₁, I ₁(u), G ₁(u) and u _LM

5) use formula (49), (46) and (25) to calculate s respectively _a, G _a(u) and u _M

6) use formula (26), (27) and (28) to calculate k respectively _G,

With

7) use respectively formula (39) and (40) calculate Γ and

8) calculate

Contrary FFT and standardization (can adjust contrast in case of necessity), promptly get restored image.

Claims (3)

1. based on the method for blindly restoring image of intermediate frequency, it is characterized in that performing step is as follows:

(1) frequency domain that contains PSF information in the image spectrum is positioned, described frequency domain localization method is:

(1.1) FFT and the standardization of calculating blurred picture;

(1.2) some of calculating FFT are crossed the frequency domain separation of former dotted line, thereby determine the scope of frequency domain;

(2) utilize frequency domain to estimate G class PSF, method is:

(2.1) according to the value of model index β in the environmental selection G class function;

(2.2) frequency domain of the identical former dotted line of mistake is carried out smoothly, and estimate the derivative of PSF in the frequency domain with the linear LSE of one-parameter;

(2.3) estimated value of fog-level parameter c in the usefulness derivative calculations G class function of PSF

Obtain the estimator of G class PSF;

(3) utilize frequency domain to optimize in the Wiener filtering noise spectral density than the value of Γ, and with Wiener filtering restored image.

2. the method for blindly restoring image based on intermediate frequency according to claim 1 is characterized in that: described step (1.2) is calculated being achieved as follows of frequency domain separation of the former dotted line of mistake of FFT:

(1.2.1) logarithm of crossing former dotted line is carried out " moving right average " smoothly, promptly the value of each point is replaced with the mean value of its several consecutive point of the right, with the straight line connection end point and calculate the point of this straight line ultimate range, promptly obtain the separation between frequency domain and the high-frequency domain then; ,

(1.2.2) " mobile left mean " carried out smoothly in the medium and low frequency territory of crossing former dotted line, promptly the value of each point is replaced with the mean value of its several consecutive point of left side, with the straight line connection end point and calculate the point of this straight line ultimate range, promptly obtain the separation between lower frequency region and the frequency domain then;

(1.2.3) zone between two separations calculating of above step is exactly this frequency domain of crossing former dotted line, when the PSF invariable rotary, then with after the frequency domain annulusization of this line as the frequency domain of image, otherwise then determine jointly by many former dotted lines of mistake.

3. the method for blindly restoring image based on intermediate frequency according to claim 1 is characterized in that: utilize frequency domain to optimize noise spectral density in the Wiener filtering in the described step (3) than the value method of Γ to be:

(3.1) utilize the middle and high frequency domain separation of the identical former dotted line of mistake to make the adaptive best Wiener filtering of Γ;

(3.2) introduce squelch parameter η, with finely tune Γ to be in harmonious proportion the contradiction of noise and image detail: η is comparatively clear near 1 o'clock noise and details, noise and details weakened simultaneously to some extent when η reduced, and therefrom chose an equilibrium value and restored with the best that realizes image.

Images