CN106157264A - Large area image uneven illumination bearing calibration based on empirical mode decomposition - Google Patents

Abstract

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition, comprise the steps: input picture, the form of detection input picture, is converted to yuv format by the image of colored rgb format, then extracts the image of Y passage in yuv format image and obtains extracting image；Convert extracting image, the image after being converted；Image after conversion is carried out decomposition based on empirical mode, obtains a series of phenogram as local frequencies and the essential mode function of dimensional properties and a surplus function；Specific surplus function and essence mode function is selected to carry out the detection of the even degree of uneven illumination；Smooth light image；Carry out reflected image conversion, obtain final correction chart picture.This method can the effective large area image under correction of complex illumination condition, real-time is good and has high subjective visual quality.

Description

Large area image uneven illumination bearing calibration based on empirical mode decomposition

Technical field

The present invention relates to image procossing and computer vision field, particularly to a kind of based on empirical mode decomposition model Large area image uneven illumination bearing calibration.

Background technology

The bearing calibration more employing homomorphic filtering of the uneven illumination of current main flow and bi-cubic interpolation method etc. are to adopting The digital picture of collection processes.Homomorphic filtering method thinks that the illumination model of image has the form of multiplication, right by doing image Transformation of variables and high-pass filtering realize, and when can not meet multiplied model for many physics illumination models, this method is imitated The most very poor.The method of bi-cubic interpolation is to think the form that illumination model has addition, by with double cubes of surface models to light Carry out regression estimates according to distribution, then the illumination deducting estimation from original image realizes, due to bi-cubic interpolation to be carried out, meter Calculation amount is very big, and interpolation precision is limited.These existing inequality illumination correction methods are mainly based upon certain physical light of priori According to model, it is impossible to accurately describe photoenvironment complicated in natural scene, and the brightness value precision estimated is relatively low, calculate complexity Spend higher, complicated scene is frequently present of robustness poor, the shortcoming that anti-noise jamming ability is weak, can not fully meet The new demand in the fields such as current aerospace shooting and satellite remote sensing.

Summary of the invention

Present invention aim to address the large area image of the complex conditions such as existing illumination acute variation and shadow occlusion The light intensity uneven distribution problem existed and existing method poor robustness, noise resisting ability difference and the problem of poor real.

To achieve these goals, the present invention provides a kind of even bearing calibration of high-precision uneven illumination to realize in complexity Under the conditions of self adaptation can remove different optical field distribution, can be applicable to large area image, and can protect with human eye subjective evaluation result Hold consistent.

Technical scheme is as follows:

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition, comprises the steps:

Input picture step, the form of detection input picture, the image of colored rgb format is converted to yuv format, then Extract the image of Y passage in yuv format image to obtain extracting image；

Image transformation step, converts extracting image, the image after being converted；

Image empirical mode decomposition step, carries out decomposition (EMD) based on empirical mode to the image after conversion, obtains one Series characterizes image local frequency and the essential mode function (IMF) of dimensional properties and a surplus function (RF), i.e.

$s\[m,n\]={\mathrm{\Σ}}_{k=0}^K{H}_k(m,n)+D,$

Wherein, the image after s [m, n] representation transformation, H_k(m n) represents the Empirical Mode collection of k essence mode function composition Closing, D represents the surplus function decomposited；

Select certain components step, select the essential mode function of k > 2 in D and Empirical Mode set to carry out uneven illumination even The detection of degree；

Smooth light image step；And

Reflected image shift step, the image obtaining described smooth light image step carries out reflected image conversion, To final correction chart picture.

Image transformation step carries out conversion to extraction image and is based on retina cerebral cortex (Retinex) algorithm, will Extracting picture breakdown and become reflecting part and luminance part, i.e. S=R × L, wherein S represents the Y channel image extracted, and R represents reflection Part, L represents luminance part,

Then S, R, L are transformed into respectively log-domain, i.e. calculate s=logS, r=logR, l=logL, wherein behalf pair The extraction image of number field, r represents the reflecting part of log-domain, and l represents the luminance part of log-domain.

In image empirical mode decomposition step, decomposition based on empirical mode includes following sub-step:

The extraction image s of log-domain is carried out EMD decomposition by row；

Use and carry out the decomposition method that EMD decomposition is identical, to log-domain with described image s that log-domain extracted by capable Extract image s and carry out EMD decomposition by row；And

The decomposition result extracting image s row and column of log-domain is carried out integrated, obtains final two-dimentional IMF and RF.

Wherein, by row, the extraction image s of log-domain is carried out EMD decomposition to comprise the steps: again

Every a line s is taken in the extraction image s of the log-domain that size is M × N_k[n] as one-dimensional input signal, if point The IMF h solved_i[n] represents, RF d_i[n] represents, initializes i=0, j=0, k=1；

Carry out screening until candidate feature signal meets two conditions of IMF, and calculate d_i=d_i-1-h_i[n]；

Repeat previous step, and set i=i+1, until stopping after meeting following condition,

$SD={\mathrm{\Σ}}_{n=0}^{N-1}\frac{{(d_i\[n\]-d_{i-1}\[n\])}^2}{d_2^i\[n\]}\∈\[a_1,a_2\],$

Wherein a₁And a₂It it is given parameter；

Repeat above-mentioned 3 steps, decompose remaining each row.

Wherein, screening process is specific as follows:

Find out s_kAll of Local modulus maxima P in [n]_maxWith local minizing point P_min, use the side of bilinear interpolation Method is according to P_maxProduce coenvelope line signal e_max[n], uses same interpolation method according to P_minProduce lower envelope line signal e_min [n]；

Calculating average tendency signal:

Calculating candidate feature signal:

J=j+1 is set, repeats above-mentioned 3 steps, until obtained c_j[n] meets two conditions of IMF.

Select the operation detecting the even degree of uneven illumination in certain components step specific as follows:

Calculate judgment criterion IR, according to equation below:

$IR={\mathrm{\σe}}^{-{\mathrm{\λ}}^2}>{\mathrm{\θ}}_1,$

Wherein, σ is picture content H_kVariance, λ is its average；θ₁For threshold value, it is usually arranged as θ₁=0.15.

If formula is satisfied, then can determine whether to there is uneven illumination even.

In smooth light image step, smooth light image includes:

(1) IMF and RF that extraction is existed uneven illumination even is smoothed, and entire image is used template size Mean filter for N × N is smoothed；

(2) use the IMF component after smoothing to replace original IMF component, then all of IMF component and RF are added Obtain the correction chart picture of log-domain.

Reflected image conversion described in reflected image shift step refers to be transformed to non-by the correction chart picture of log-domain The correction chart picture of log-domain；For coloured image, also it is transformed to RGB territory from YUV territory.

Beneficial effects of the present invention is as follows:

The present invention relates to a kind of method utilizing empirical mode decomposition to carry out image uneven illumination correction.The present invention passes through Image is carried out multiple dimensioned and frequency and combines decomposition, utilize the uneven light intensity of design and shade existence rules to judge, Obtain the accurate estimation of light distribution in image, after smoothing processing, from original image, finally deduct plot of light intensity just can obtain Image after correction.This method can the effective large area image under correction of complex illumination condition, can calculate in real time, and And have high subjective visual quality.

Accompanying drawing explanation

Fig. 1 is the flow process of present invention large area image based on empirical mode decomposition model uneven illumination bearing calibration Figure；

Fig. 2 is to extracting the algorithm principle figure that image converts in the embodiment of the present invention；

When Fig. 3 is in the embodiment of the present invention to convert extraction image, the processing procedure of the Retinex algorithm of foundation is shown It is intended to.

Fig. 4 is the applicating flow chart of EMD method in the embodiment of the present invention；

Fig. 5 is that wherein (a) is the image before correcting according to the effect contrast figure of image before and after the method correction of the present invention, B () is the image after correction.

Fig. 6 is according to Y-PSNR (PSNR) figure of image before and after the method correction of the present invention, before wherein (a) is for correction Image PSNR figure, (b) be correction after image PSNR scheme.

Detailed description of the invention

For making the object, technical solutions and advantages of the present invention clearer, below in conjunction with concrete example, and with reference to attached Figure, the invention will be further described.

Hardware environment for implementing is: computer CPU is Intel Core I5, inside saves as 4GB.Software environment is: Matlab R2010b and Windows7.

The large area image uneven illumination bearing calibration based on empirical mode decomposition model of the present invention utilizes EMD method Image is carried out the multi-scale Representation of data adaptive, obtain the IMF of a series of expression image local frequency and dimensional properties with And a RF, choose qualified IMF and the RF estimation as image illumination component according to the judgment criterion of design, then to this After a little luminance component use smoothing operator to process, remove from original input picture, it is possible to finally give the school of uniform illumination Positive image.As it is shown in figure 1, in short, 6 steps are altogether included, it may be assumed that input picture；Image converts；Image empirical mode decomposition； Select certain components；Smooth light image；Reflected image converts.By above step, the image after being corrected.

The present invention is embodied as follows:

In input picture step S1, carry out input image data preparation: detect the form of input picture, if colored The image of rgb format is transformed into yuv format, and the image extracting Y passage obtains extracting image.

Rgb format is converted to yuv format to be carried out according to following equation:

$\left[\begin{array}{c}Y\\ U\\ V\end{array}\right]=\left[\begin{array}{ccc}0.299& 0.587& 0.114\\ -0.169& -0.331& 0.5\\ 0.5& -0.419& -0.081\end{array}\right]\left[\begin{array}{c}R\\ G\\ B\end{array}\right]+\left[\begin{array}{c}0\\ 128\\ 128\end{array}\right],$

In image transformation step S2, convert extracting image, specific as follows:

According to Retinex algorithm, decomposing extraction image has, S=R × L, and wherein S represents the Y channel image extracted, and R represents Reflecting part, L represents luminance part,

Then S, R, L are transformed into respectively log-domain, i.e. calculate s=logS, r=logR, l=logL, wherein behalf pair The extraction image of number field, r represents the reflecting part of log-domain, and l represents the luminance part of log-domain.

Image S carrying out conversion theoretical according to Retinex in the present embodiment, this algorithm is thought: the image S that a width is given (x, y) can be decomposed into two different images: reflected image R (x, y) and luminance picture (also referred to as incident image) L (x, Y), i.e. S (x, y)=R (and x, y) L (x, y), its principle is: image can be regarded as being made up of incident image and reflected image, Incident illumination is radiated on reflection object, enters human eye by the reflection light that reflects to form of reflection object, it is simply that the mankind are seen Image, as shown in Figure 2.

Wherein, (x y) directly determines the dynamic range that in image, pixel can reach to L；(x y) represents the anti-of object to R Penetrate the inherent attribute of character image, i.e. image；(x y) represents the reflected light image that human eye can receive to S.Retinex is theoretical Basic thought be exactly in original image, by removing or reduce the impact of incident image someway, thus as best one can Retain the reflecting attribute image of material essence.The general processing procedure schematic diagram of Retinex algorithm is as shown in Figure 3.

The formula of Retinex algorithm is:

S (x, y)=R(x, y) L(x, y),

$r(x,y)=\log R(x,y)=log\frac{S(x,y)}{L(x,y)},$

R(x, y)=logR (x, y)-log [F (and x, y) * S (x, y)],

Here, r (x, y) for output image, * is convolution symbol, F (x, y) centered by around function, can be expressed as:

$F(x,y)={\mathrm{\λe}}^{\frac{-(x^2+y^2)}{C^2}},$

Wherein, C is Gauss around yardstick, and λ is a yardstick, and its value meets following condition:

∫ ∫ F (x, y) dxdy=1.

Retinex algorithm to realize flow process as follows:

(1) read in artwork, if gray-scale map, then the gray value of each for image pixel be converted to floating number by integer value, And be transformed into log-domain, if input be coloured picture, each color classification of image is processed respectively, by the pixel value of each component by Integer value is converted to floating number, and is transformed in log-domain；

(2) input yardstick C；Under discrete conditions, integration is converted into summation；Further determine that the value of parameter lambda；

(3) according to formula above, be calculated r (x, y)；If coloured picture, the most each passage all has a r_i(x, y)；

(4) by r (x, y) is transformed into real number field from log-domain, obtain export image R (x, y)；

(5) to R (x, y) linear stretch with corresponding form output display.

Image empirical mode decomposition is carried out, its sub-process schematic diagram such as Fig. 4 institute in image empirical mode decomposition step S3 Show, specific as follows:

In step S31, the extraction image of the log-domain that size is M × N is taken every a line s_k[n] is as one-dimensional input Signal, if the IMF h decomposited_i[n] represents, RF d_i[n] represents, initializes i=0, j=0, k=1.

In step s 32, s is found out_kAll of Local modulus maxima P in [n]_maxWith local minizing point P_min, use double The method of linear interpolation is according to P_maxProduce coenvelope line signal e_max[n], uses same interpolation method according to P_minProduce lower bag Winding thread signal e_min[n]。

In step S33, calculating average tendency signal:

In step S34, calculating candidate feature signal:

In step s 35, j=j+1 is set, repeats step 4-6, until obtained c_j[n] meets two bars of IMF Part.

In step S36, h is set_i[n]=c_j[n], calculating RF has d_i=d_i-1-h_i[n]

In step S37, repeat step 4-8, and set i=i+1, until meeting following condition stopping SD:

$SD={\mathrm{\Σ}}_{n=0}^{N-1}\frac{(d_i\[n\])-d_{i-1}{\[n\]}^2}{d_2^i\[n\]}\∈\[0.2,0.3\].$

Step S31-step S37 is the EMD process of a standard, and EMD decomposition method is based on it is assumed hereinafter that condition: (1) Data at least two extreme values, a maximum and a minima；(2) between the local temporal characteristic of data is by extreme point Time scale uniquely determines；(3) if data do not have extreme point but have flex point, then can be by Numeric differential one or many Try to achieve extreme value, obtain decomposition result by integration the most again.The essence of this method is the characteristic time scale by data Obtain intrinsic fluctuation model, then decomposition data.This catabolic process is referred to as " screening " process.

Two conditions of IMF are: (1) function is in whole time range, and the number of Local Extremum and zero crossing is necessary Equal, or at most difference one；(2) point, envelope (the coenvelope line e of local maximum at any time_max[n]) and Envelope (the lower envelope line e of little value_min[n]) averagely it is necessary for zero.If the c that screening obtains_j[n] is unsatisfactory for the two condition, then Needs proceed screening.

The implementation process of the EMD decomposition method of the present embodiment sees Fig. 4.

Use identical step S31-step S37, decompose remaining each row.The extraction image of log-domain is pressed by same mode Row carry out EMD decomposition, finally give the result that M × N number of one-dimensional signal decomposes.

The decomposition result of row and column is carried out integrated, final two-dimentional IMF and RF can be obtained, the EMD table of two dimensional image It is shown as:

$s\[m,n\]={\mathrm{\Σ}}_{k=0}^K{H}_k(m,n)+D.$

By step S3, complete the EMD process of two dimensional image.This process is based on an assumption that (1) 2-D data plane There is no extreme point including at least a maximum point and a minimum point or whole two dimensional surface but carrying out single order or several rank One maximum point and a minimum point can occur after derivative operation；(2) chi of the spacing of characteristic dimension extreme point Degree definition.

In the present embodiment, by the EMD process of two dimensional image, original function is broken down into the Empirical Mode set of k IMF composition (H_k(m, n)) and remainder D (i.e. RF).Wherein, H_k(m is n) the different scale figure layer obtained after Scale separation, with k's Increasing, yardstick increases；D is final trend term.

In selecting certain components step S4, choose certain components and carry out the even degree detecting of uneven illumination.Due in image The details such as texture be primarily present in high-frequency IMF component, therefore to k > the two-dimentional H of 2_kThe even journey of uneven illumination is carried out with D Degree and the detection of shade, it is judged that criterion IR is:

$IR={\mathrm{\σe}}^{-{\mathrm{\λ}}^2}>{\mathrm{\θ}}_1,$

Wherein, σ is picture content H_kVariance, λ is its average.Threshold value is set to θ₁=0.15.If above formula meets, Then can determine whether to there is uneven illumination even.

In smooth light image step S5, light image is smoothed, including:

(1) IMF and D that extraction is existed uneven illumination even is smoothed, and entire image is used template size Be 30 × 30 mean filter be smoothed.

(2) use the IMF component after smoothing to replace original IMF component, then all of IMF component and D are added Obtain the correction chart picture of log-domain.

In image reflection transformation step S6, finally the correction chart picture in log-domain is transformed to the correction in non-logarithmic territory Image, if coloured image, then transform to RGB territory from YUV territory.

Through step S1-S6, obtain correction chart picture.

Fig. 5 show the image effect comparison diagram before and after using the method for the present invention to carry out illumination correction, and wherein (a) is school Before just, after (b) is for correction, it can be seen that before and after using the method for the present invention to carry out illumination correction, the light distribution of image Uniformity coefficient is greatly promoted.

Fig. 6 show the image PSNR figure before and after using the method for the present invention to carry out illumination correction, before wherein (a) is for correction , after (b) is for correction, it can be seen that before and after using the method for the present invention to carry out illumination correction, the PSNR value of image drops significantly Low, it was demonstrated that its noise resisting ability of the image after corrected is also greatly enhanced.

Claims (10)

1. a large area image uneven illumination bearing calibration based on empirical mode decomposition, it is characterised in that include as follows Step:

Input picture step, the form of detection input picture, the image of colored rgb format is converted to yuv format, then extracts In yuv format image, the image of Y passage obtains extracting image；

Image transformation step, converts described extraction image, the image after being converted；

Image empirical mode decomposition step, carries out decomposition based on empirical mode, obtains a series of phenogram the image after described conversion As local frequencies and the essential mode function of dimensional properties and a surplus function, i.e. Wherein, the image after s [m, n] representation transformation, H_k(m, n) represents the Empirical Mode set of k essence mode function composition, and D represents The surplus function decomposited；

Select certain components step, select the essential mode function of k > 2 in surplus function and Empirical Mode set to carry out illumination not The detection of uniformity coefficient；

Smooth light image step；And

Reflected image shift step, the image obtaining described smooth light image step carries out reflected image conversion, obtains Whole correction chart picture.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 1, It is characterized in that, described in described image transformation step, extraction image is carried out conversion and refers to calculate according to retina cerebral cortex Method, becomes reflecting part and luminance part, i.e. S=R × L by extraction picture breakdown, and wherein S represents the Y channel image extracted, R generation Table reflecting part, L represents luminance part,

Then S, R, L are transformed into respectively log-domain, i.e. calculate s=logS, r=logR, l=logL, wherein behalf log-domain Extraction image, r represents the reflecting part of log-domain, and l represents the luminance part of log-domain.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 1, It is characterized in that, decomposition based on empirical mode described in described image empirical mode decomposition step includes following sub-step:

The extraction image of log-domain is carried out empirical mode decomposition by row；

Use, with described, the image that extracts of log-domain is carried out decomposition method identical in empirical mode decomposition, to log-domain by row Extraction image carry out empirical mode decomposition by row；And

The decomposition result extracting image row and column of log-domain is carried out integrated, obtain final two dimension essential mode function and One surplus function.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 3, It is characterized in that, described the extraction image of log-domain is carried out empirical mode decomposition is comprised the steps: by row

Every a line s is taken in the extraction image s of the log-domain that size is M × N_k[n] as one-dimensional input signal, if decompositing Essential mode function h_i[n] represents, surplus function d_i[n] represents, initializes i=0, j=0, k=1；

Carry out screening until candidate feature signal meets two conditions of essence mode function, and calculate d_i=d_i-1-h_i[n]；

Repeat previous step, and set i=i+1, until stopping after meeting following condition,

$SD={\mathrm{\Σ}}_{n=0}^{N-1}\frac{{(d_i\[n\]-d_{i-1}\[n\])}^2}{d_2^i\[n\]}\∈\[{\mathrm{\α}}_1,{\mathrm{\α}}_2\],$

Wherein a₁And a₂It it is given parameter；

Repeat above-mentioned 3 steps, decompose remaining each row.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 4, It is characterized in that, described screening process is as follows:

Find out described s_kAll of Local modulus maxima P in [n]_maxWith local minizing point P_min, the method for use bilinear interpolation According to Local modulus maxima P_maxProduce coenvelope line signal e_max[n], uses same interpolation method according to local minizing point P_minProduce lower envelope line signal e_min[n]；

Calculate average tendency signal

Calculate candidate feature signalAnd

J=j+1 is set, repeats above-mentioned 3 steps, until obtained candidate feature signal meets two of essence mode function Condition.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 1, It is characterized in that, the operation detecting the even degree of uneven illumination described in described selection certain components step is specific as follows:

Calculate judgment criterion IR, according to equation below:

$IR={\mathrm{\σe}}^{-{\mathrm{\λ}}^2}>{\mathrm{\θ}}_1,$

Wherein, σ is picture content H_kVariance, λ is its average, θ₁For threshold value；

If meeting described formula, then it is judged as there is uneven illumination even.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 6, It is characterized in that, described threshold value is 0.15.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 1, It is characterized in that, the specific practice of described smooth light image step is:

Be there is the even essential mode function of uneven illumination in extraction and surplus function is smoothed, entire image is used Template size is that the mean filter of N × N is smoothed；And

The essential mode function after smoothing is used to replace the essential mode function before smoothing, then by all of essence pattern letter Number and surplus function are added the correction chart picture obtaining log-domain.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 1, It is characterized in that, the reflected image conversion described in described reflected image shift step is to be transformed to by the correction chart picture of log-domain The correction chart picture in non-logarithmic territory.

A kind of large area image uneven illumination bearing calibration based on empirical mode decomposition the most according to claim 9, It is characterized in that, the conversion of described reflected image also includes, for coloured image, from YUV territory, it being transformed to RGB territory.