CN101364302A - Clearness processing method for defocus blurred image - Google Patents

Abstract

The invention discloses a defocus blurred image sharpness processing method, which is implemented according to the following steps: firstly, the mean value and the variance of the edge width is calculated according to gradient information of an image; an initial parameter of a fuzzy model is obtained according to the statistical data; secondly, the blurred image is divided into concentric circles around the center of the image and is broken down into sub-images, wherein, the number of the concentric circles is k and the number of the sub-images is k+1; the corresponding blurred initial semidiameters are distributed to the sub-images; the optimal blurred semidiameter is found through an iterative mode, and the sharpness processing of all sub-images is performed by adopting a frequency domain inverse-filtering manner; finally, all the sub-images are added to synthesize the whole sharpness image. The defocus blurred image sharpness processing method overcomes the limitations in the existing circular disc function modeling recovery method based on the fixed semidiameter, and serves the purpose for restoring the sharp image.

Description

A kind of clearness processing method of defocus blurred image

Technical field

The invention belongs to digital picture recovery technology field, relate to a kind of clearness processing method of defocus blurred image.

Background technology

When taking pictures, image blurring situation appears sometimes, and the reason that causes this situation mainly is that subject is not on the focal plane of imaging system, is referred to as defocusing blurring.For some photos that comprises important information, can't discern because of defocusing blurring, if having no chance to take once more, it is too high perhaps to take the cost that will pay once more, just can adopt the recovery sharpening technology of defocusing blurring that it is carried out the processing of sharpening.The recovery sharpening technology of existing defocusing blurring is based on radii fixus disk function modelling restoration methods mostly, and the recovery capability that these class methods are blured the optical defocus that produces in the reality is extremely limited.

Summary of the invention

The object of the present invention is to provide a kind of clearness processing method of defocus blurred image, overcome prior art, the image of defocusing blurring can be reverted to distinct image the limited problem of the recovery capability of blurred picture.

The technical solution adopted in the present invention is, a kind of clearness processing method of defocus blurred image, and this method is implemented according to following steps,

Step 1: pending defocus blurred image is set up mathematical model

Utilize a disk function to carry out the modeling of blurred picture, that is:

$h(x,y)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\π}{R_b}^2}& x^2+y^2\≤{R_b}^2\\ 0& x^2+y^2>{R_b}^2\end{array}\right.---\left(1\right)$

Wherein, R _bBe blur radius, π is a circular constant, and (x y) is certain pixel on the image;

Step 2: choose initial blur radius

Whole figure carries out searching of blur radius to pending defocusing blurring, at first taken out-of-focus image is carried out sharpening, generates gradient map E _sWith directional diagram E _oTo the gradient map E that obtains _sAnd gradient direction figure E _o, obtain partial gradient amplitude maximum value point set E with non-maximum value inhibition technology _nTo E _nUnder given threshold value, carry out binaryzation, obtain profile diagram E _e, ask for the edge width again, the width value that obtains is calculated average M and variance D, choose the initial blur radius R of M/2 for whole figure _b

Step 3: the division of subimage

With resulting variance D substitution equation k=int[D-1 in the step 2], obtain the parameter k that subgraph is divided number, be round dot again with the picture centre, draw k concentric circles, its radius is R _k=R ₁+ (k-1) Δ R, wherein R ₁Be the innermost circle radius, Δ R is an increment, is a border circular areas, a k-1 circle ring area, residual image zone k+1 number of sub images altogether with picture breakdown, is designated as g _S1, g _S2... g _{S (k+1)}

Step 4: the sharpening of subimage is handled

With the k that obtains in variance that obtains in the step 2 and the step 3, the step-length that obtains disc radius is D/k, obtains the initial blur radius under the i number of sub images, according to this blur radius, obtains the disk function h of subimage according to equation (1) _i(x, y), with subimage g _Si(x y) carries out the frequency domain conversion, again disk function corresponding under each subimage is carried out frequency domain transform, and the power spectrum S of calculating noise and original image _Nn(u, v) and S _Ff(u, v); Carry out clear picture based on liftering then and recover, change blur radius, promptly calculate R _i(0) and R _i(0) ± Δ R _iObtain three with this and recover subimage, calculate these three the Sobel sharpen detail energygrams that recover image, and calculate these three differences of recovering the image variances, carry out iteration according to the little direction of difference then, on the opportunity up to searching out variance difference minimum, stop iteration, the image of the ceiling capacity average in two energygrams when selecting iteration stopping is a net result, is made as f _i* (x, y);

Step 5, acquisition view picture sharpening image

With f _i* (x y) carries out addition,

Promptly get the sharpening image $f(x,y)=\underset{i=1}{\overset{k+1}{\mathrm{\Σ}}}{f}_i*(x,y).$

The clearness processing method of defocus blurred image of the present invention; overcome existing limitation based on radii fixus disk function modelling restoration methods; utilize the basic image-forming principle of defocus blurred image, build the degeneration fuzzy model disk function that becomes radius, reach the purpose that recovers picture rich in detail.

Description of drawings

Fig. 1 is for defocusing principle schematic, and wherein a is the ornaments relation of camera and photographed, and b is for defocusing principal diagram;

Fig. 2 is the subimage decomposing schematic representation in the inventive method, and wherein a is concentrically ringed division, and b is the regional extent of k subgraph.

Embodiment

The present invention is described in detail below in conjunction with the drawings and specific embodiments.

The present invention proposes a kind of restoration methods that becomes under the radius disk function, it is basic image-forming principle according to defocus blurred image, make up the degeneration fuzzy model disk function of different radii, blurred picture is decomposed into a plurality of subimages, each subimage is searched best blur radius by alternative manner respectively, and the method for employing frequency domain inverse filtering, the defocusing blurring subimage is carried out liftering, thereby reach the purpose that recovers picture rich in detail.

The clearness processing method of defocus blurred image of the present invention, implement according to the following steps:

Step 1:, set up mathematical model according to the degeneration principle of defocus blurred image

Fig. 1 a is the relative position synoptic diagram of digital image-forming equipment and photographic subjects thing, and this is a kind of desirable shooting state.Shown in Fig. 1 b, when the imaging plate (CCD plane) that does not project camera inside when picture point is gone up, just can produce defocusing blurring.The light cone path that on behalf of the diverse location by object, four dotted lines among Fig. 1 b send in twos, picture point is that the summit of light cone is not dropped on the camera imaging plate, what therefore accept on the imaging plate will be the xsect of light cone, i.e. hot spot, image is formed by stacking by these hot spots, so form blur effect.The shape that the camera imaging plate is accepted the sub-circular hot spot and radius can be because different with the relative position (being the center of image) of optical center and change, its Changing Pattern is for approaching picture centre, hot spot just approaches circle more, and the radius of hot spot is just more little; Away from picture centre, spot radius increases progressively more.

Therefore, defocus can be approximated to be by some expansion for an equally distributed circular light spot, and radius so just can be reduced to it disk function and carry out the modeling of blurred picture, that is: along with increasing progressively away from picture centre

$h(x,y)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\π}{R_b}^2}& x^2+y^2\≤{R_b}^2\\ 0& x^2+y^2>{R_b}^2\end{array}\right.---\left(1\right)$

Wherein, R _bBe blur radius, π is a circular constant, and (x y) is certain pixel on the image.

Step 2: the choosing of initial blur radius

The mathematical model that has obtained defocusing blurring by step 1 as can be known, blur radius R _bBe the parameter of determining fog-level, therefore, earlier whole figure carried out searching of blur radius according to following method, the concrete steps of its acquisition are as follows:

2.1, the image gradient information extraction

Taken out-of-focus image is carried out sharpening, generate gradient map E _sWith gradient direction figure E _o, sharpening adopts the Sobel operator, and the calculation equation of gradient is:

$S(i,j)=\sqrt{d_x^2+d_y^2}---\left(2\right)$

Wherein,

G represents out-of-focus image.

Gradient map E then _sBe E _s=[S (i, j)].

The gradient direction angle, ask for by following formula:

${\mathrm{\θ}}_M(i,j)={\mathrm{tg}}^{-1}\left(\frac{d_y}{d_x}\right)---\left(3\right)$

To the θ that asks for _MBy 90 °～67.5 °, 67.5 °～22.5 °, 22.5 °～-22.5 ° ,-22.5 °～-67.5 ° ,-67.5 °～-90 ° five scopes are divided into 90 °, and 45 °, 0 ° ,-45 ° ,-90 ° of five directions constitute gradient direction figure E _o

2.2, obtain partial gradient amplitude maximum value point set E with non-maximum value inhibition technology _n

The gradient map E that step 2.1 is obtained _s, and gradient direction figure E _o, to each pixel g (x, y), according to gradient direction figure E _oIn pointed direction d _k, along two adjacent pixels of positive and negative biplane examinations.If E _s(x y) greater than the gradient intensity of two neighbors, so just makes E _n(x, y)=E _s(x y), otherwise makes E _n(x, y)=0.The matrix E that obtains at last _n(x y) is partial gradient amplitude maximum value point set.

2.3, the asking for of edge width

To local gradient amplitude maximum value point set E _n, under given threshold value, carry out binaryzation, obtain profile diagram E _eAfterwards, to E _nIn any 1 P, according to gradient direction figure E _oIndicated direction is along the positive and negative direction search E of this point _eIn the point.If along positive and negative direction search to E _eIn first point be respectively P ₁, P ₂If, P ₁And P ₂All satisfy: ‖ P-P ₁‖≤d and ‖ P-P ₂‖≤d (wherein d is given threshold value, i.e. Sou Suo longest distance) is then with P ₁And P ₂Between distance as the respective edges width, otherwise skip this point.Average M, variance D are tried to achieve in the width value calculating that obtains.Choose the initial blur radius R of M/2 for whole figure _b

Step 3: the division of subimage

When considering imaging, the fog-level that defocuses is different on entire image, so the present invention is decomposed into a plurality of subimages with original image.

Decomposition method is: shown in Fig. 2 a, be round dot with the picture centre, draw k concentric circles, its radius is R _k=R ₁+ (k-1) Δ R.R wherein ₁Be the innermost circle radius, Δ R is an increment.So just can be a border circular areas, a k-1 circle ring area, residual image zone k+1 number of sub images altogether with picture breakdown, be designated as g _S1, g _S2... g _{S (k+1)}, the dash area of the regional extent of k subgraph shown in Fig. 2 b.

Step 4: the sharpening of subimage is handled

The step-length that can be got disc radius by the variance that obtains in the step 2 is D/k, and the initial blur radius that defines under the i number of sub images can be expressed as so:

R _bi(0)＝R _b+(i-1)·D/k (i＝1，2，3...，k+1)。

4.1, according to this blur radius, can obtain the disk function h of subimage according to equation (1) _i(x, y)

4.2, with blurred picture subimage g _Si(x y) utilizes discrete two-dimensional Fourier transform (DFT) to finish the frequency domain conversion respectively, that is:

$G_i(u,v)=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{g}_{\mathrm{si}}(x,y)\exp [-\frac{2\mathrm{\πj}}{N}(\mathrm{ux}+\mathrm{vy})]---\left(4\right)$

Wherein, G _i(u, v), i=1,2,3...k+1 is the frequency domain transform of subimage.

4.3, the disk function that each subimage is down corresponding carries out frequency domain transform.That is:

$H_i(u,v)=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{h}_i(x,y)\exp [-\frac{2\mathrm{\πj}}{N}(\mathrm{ux}+\mathrm{vy})]---\left(5\right)$

Wherein, h _i(u v) is the corresponding disk function of subimage, H _i(u, v) system function behind the frequency domain transform.

4.4, the power spectrum S of calculating noise and original image _Nn(u, v) and S _Ff(u, v)

Directly calculate near the local variance of the collection of pixels each pixel, choose maximal value in the local variance, on image, look for a flat site simultaneously as the variance of image from blurred picture, with its local variance as noise variance.But often manually be not easy to find flat site, so can utilize the local variance of following formula computed image, the image boundary variance is not taken into account.With the estimation of the ratio of the maximal value of local variance and minimum value, that is: as signal noise ratio (snr) of image

${{\mathrm{\σ}}^2}_{\mathrm{yL}}(i,j)=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}{[y(i+k,j+l)-{\mathrm{\μ}}_y(i,j)]}^2---\left(6\right)$

Wherein, μ _yBe local mean value, be calculated as follows:

${\mathrm{\μ}}_y=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}y(i+k,j+l)---\left(7\right)$

The size that variance is calculated the window that uses is P=Q=2 (i.e. 5 * 5 windows).

4.5, recover based on the clear picture of liftering

According to the rapid subimage function H of previous step _i(u v) calculates each complex conjugate function The present invention selects for use S filter to carry out liftering and handles, thereby realizes sharpening, and the frequency domain presentation of each subimage sharpening image is designated as F _i(u, v), that is:

$F_i(u,v)=\frac{{H^{*}}_i(u,v)G_i(u,v)}{{\left|H_i(u,v)\right|}^2+S_{\mathrm{nn}}(u,v)/S_{\mathrm{ff}}(u,v)}---\left(8\right)$

To each subimage F under the frequency domain _i(u v) carries out two-dimensional discrete Fu Shi inverse transformation, recovers the original image f of each subimage _i(x, y).

$f_i(x,y)=\frac{1}{N^2}\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{F}_i(u,v)\exp \left[\frac{2\mathrm{\πj}}{N}(\mathrm{ux}+\mathrm{vy})\right]---\left(9\right)$

4.6, change blur radius, promptly calculate R _Bi(0) and R _Bi(0) ± Δ R _i(Δ R _iSuggestion is 1 pixel), according to step 4.1～step 4.5 acquisition recovery subimage separately.

4.7, obtain three Sobel sharpen detail energygrams that recover images of calculation procedure 4.6, the acquisition methods of energygram adopts equation (2).

4.8, obtain three differences of recovering the image variances of calculation procedure 4.6, carry out iteration according to the little direction of difference then, on the opportunity up to searching out variance difference minimum, stop iteration, the image of the ceiling capacity average in two energygrams when selecting iteration stopping is a net result, is made as f _i* (x, y).So far, the sharpening of finishing all subimages recovers to handle.

Step 5: obtain view picture sharpening image

With f _i* (x y) carries out addition,

$f(x,y)=\underset{i=1}{\overset{k+1}{\mathrm{\Σ}}}{f}_i*(x,y),$ Promptly obtain the sharpening image.

Claims (4)

1. the clearness processing method of a defocus blurred image is characterized in that, this method is implemented according to following steps,

Step 1: pending defocus blurred image is set up mathematical model

Utilize a disk function to carry out the modeling of blurred picture, that is:

$h(x,y)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\π}{R_b}^2}& x^2+y^2\≤{R_b}^2\\ 0& x^2+y^2>{R_b}^2\end{array}\right.---\left(1\right)$

Wherein, R _bBe blur radius, π is a circular constant, and (x y) is certain pixel on the image;

Step 2; Choose initial blur radius

Whole figure carries out searching of blur radius to pending defocusing blurring, at first taken out-of-focus image is carried out sharpening, generates gradient map E _sWith directional diagram E _oTo the gradient map E that obtains _sAnd gradient direction figure E _o, obtain partial gradient amplitude maximum value point set E with non-maximum value inhibition technology _nTo E _nUnder given threshold value, carry out binaryzation, obtain profile diagram E _e, ask for the edge width again, the width value that obtains is calculated average M and variance D, choose the initial blur radius R of M/2 for whole figure _b

Step 3: the division of subimage

With resulting variance D substitution equation k=int[D-1 in the step 2], obtain the parameter k that subgraph is divided number, be round dot again with the picture centre, draw k concentric circles, its radius is R _k=R ₁+ (k-1) Δ R, wherein R ₁Be the innermost circle radius, Δ R is an increment, is a border circular areas, a k-1 circle ring area, residual image zone k+1 number of sub images altogether with picture breakdown, is designated as g _S1, g _S2... g _s(k+1);

Step 4: the sharpening of subimage is handled

With the k that obtains in variance that obtains in the step 2 and the step 3, the step-length that obtains disc radius is D/k, obtains the initial blur radius under the i number of sub images, according to this blur radius, obtains the disk function h of subimage according to equation (1) _i(x, y), with subimage g _Si(x y) carries out the frequency domain conversion, again disk function corresponding under each subimage is carried out frequency domain transform, and the power spectrum S of calculating noise and original image _Nn(u, v) and S _Ff(u, v); Carry out clear picture based on liftering then and recover, change blur radius, promptly calculate R _i(0) and R _i(0) ± Δ R _iObtain three with this and recover subimage, calculate these three the Sobel sharpen detail energygrams that recover image, and calculate these three differences of recovering the image variances, carry out iteration according to the little direction of difference then, on the opportunity up to searching out variance difference minimum, stop iteration, the image of the ceiling capacity average in two energygrams when selecting iteration stopping is a net result, is made as f _i ^*(x, y);

Step 5, acquisition view picture sharpening image

With f _i ^*(x y) carries out addition,

Promptly get the sharpening image $f(x,y)=\underset{i=1}{\overset{k+1}{\mathrm{\Σ}}}{f}_i*(x,y).$

2. according to the described clearness processing method of claim 1, it is characterized in that the described initial blur radius of choosing is taked following concrete steps:

2.1, the image gradient information extraction

Taken out-of-focus image is carried out sharpening, generate gradient map E _sWith gradient direction figure E _o, sharpening adopts the Sobel operator, and the calculation equation of gradient is:

$S(i,j)=\sqrt{d_x^2+d_y^2}---\left(2\right)$

Wherein,

,

, g represents out-of-focus image, then gradient map E _sBe E _s=[S (i, j)],

The gradient direction angle, ask for by following formula:

${\mathrm{\θ}}_M(i,j)={\mathrm{tg}}^{-1}\left(\frac{d_y}{d_x}\right)---\left(3\right)$

To ask for θ _MBy 90 °～67.5 °, 67.5 °～22.5 °, 22.5 °～-22.5 ° ,-22.5 °～-67.5 ° ,-67.5 °～-90 ° five scopes are divided into 90 °, and 45 °, 0 ° ,-45 ° ,-90 ° of five directions constitute gradient direction figure E ₀

2.2, obtain partial gradient amplitude maximum value point set En with non-maximum value inhibition technology

The gradient map E that step 2.1 is obtained _s, and gradient direction figure E _o, to each pixel g (x, y), according to gradient direction figure E _oIn pointed direction d _k, along two adjacent pixels of positive and negative biplane examinations, if E _s(x y) greater than the gradient intensity of two neighbors, so just makes E _n(x, y)=E _s(x y), otherwise makes E _n(x, y)=0, the matrix E that obtains at last _n(x y) is partial gradient amplitude maximum value point set;

2.3, the asking for of edge width

To local gradient amplitude maximum value point set E _n, under given threshold value, carry out binaryzation, obtain profile diagram E _e, then to E _nIn any 1 P, according to gradient direction figure E _oIndicated direction is along the positive and negative direction search E of this point _eIn point, establish along positive and negative direction search to E _eIn first point be respectively P ₁, P ₂If, P ₁And P ₂All satisfy: || P-P ₁||≤d and || P-P ₂||≤d, wherein d is given threshold value, i.e. Sou Suo longest distance is then with P ₁And P ₂Between distance as the respective edges width, otherwise skip this point,

Average M, variance D are tried to achieve in the width value calculating that obtains, choose the initial blur radius R of M/2 for whole figure _b

3. according to the described disposal route of claim 1, it is characterized in that the sharpening of subimage is handled according to following concrete steps and implemented in the described step 4,

By the k that obtains in variance that obtains in claim 1 step 2 and the step 3, the step-length that obtains disc radius is D/k, and the initial blur radius that obtains under the i number of sub images is:

R _bi(0)＝R _b+(i-1)·D/k (i＝1，2，3...，k+1)

4.1, according to this blur radius, obtain the disk function h of subimage according to equation (1) _i(x, y);

4.2, with the subimage g of blurred picture _Si(x y) utilizes the discrete two-dimensional Fourier transform to carry out the frequency domain conversion respectively, that is:

$G_i(u,v)=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{g}_{\mathrm{si}}(x,y)\exp [-\frac{2\mathrm{\πj}}{N}(\mathrm{ux}+\mathrm{vy})]---\left(4\right)$

Wherein, G _i(u, v), i=1,2,3...k+1 is the frequency domain transform of subimage;

4.3, the disk function that each subimage is down corresponding carries out frequency domain transform, that is:

$H_i(u,v)=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{h}_i(x,y)\exp [-\frac{2\mathrm{\πj}}{N}(\mathrm{ux}+\mathrm{vy})]---\left(5\right)$

Wherein, h _i(u v) is the corresponding disk function of subimage, H _i(u, v) system function behind the frequency domain transform;

4.4, the power spectrum S of calculating noise and original image _Nn(u, v) and S _Ff(u, v);

Near the local variance of direct collection of pixels each pixel of blurred picture calculating, choose maximal value in the local variance as the variance of image, on image, look for simultaneously a flat site, with its local variance as noise variance, utilize the local variance of equation (6) computed image, with the estimation of the ratio of the maximal value of local variance and minimum value, that is: as signal noise ratio (snr) of image

${{\mathrm{\σ}}^2}_{\mathrm{yL}}(i,j)=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}{[y(i+k,j+l)-{\mathrm{\μ}}_y(i,j)]}^2---\left(6\right)$

μ in the equation (6) _yBe local mean value, be calculated as follows:

${\mathrm{\μ}}_y=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}y(i+k,j+l)---\left(7\right)$

The size that variance is calculated the window that uses is P=Q=2;

4.5, recover based on the clear picture of liftering

Subimage function H according to step 4.3 _i(u v) calculates each complex conjugate function Carry out liftering with wave filter and handle, thereby realize sharpening, the frequency domain presentation of each subimage sharpening image is designated as F _i(u, v), that is:

$F_i(u,v)=\frac{{H^{*}}_i(u,v)G_i(u,v)}{{\left|H_i(u,v)\right|}^2+S_{\mathrm{nn}}(u,v)/S_{\mathrm{ff}}(u,v)}---\left(8\right)$

To each subimage F under the frequency domain _i(u v) carries out two-dimensional discrete Fu Shi inverse transformation, recovers the original image f of each subimage _i(x, y),

$f_i(x,y)=\frac{1}{N^2}\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{F}_i(u,v)\exp \left[\frac{2\mathrm{\πj}}{N}(\mathrm{ux}+\mathrm{vy})\right]---\left(9\right)$

4.6, change blur radius, promptly calculate R _i(0) and R _i(0) ± Δ R _i, obtain three according to step 4.1 to step 4.5 and recover subimage;

4.7, obtain three Sobel sharpen detail energygrams that recover images of calculation procedure 4.6, the acquisition methods of energygram adopts equation (2);

4.8, obtain three differences of recovering the image variances of calculation procedure 4.6, afterwards, carry out iteration according to the direction that difference is little, up to the opportunity that searches out variance difference minimum, stop iteration, the image of the ceiling capacity average in two energygrams when selecting iteration stopping is a net result, is made as f _i ^*(x, y).

Images