CN102609904A - Bivariate nonlocal average filtering de-noising method for X-ray image - Google Patents

Abstract

The invention provides a bivariate nonlocal average filtering de-noising method for an X-ray image. The method is characterized by comprising the following steps: 1) a selecting method of a fuzzy de-noising window; and 2) a bivariate fuzzy adaptive nonlocal average filtering algorithm. The method has the beneficial effects that in order to preferably remove the influence caused by the unknown quantum noise existing in an industrial X-ray scan image, the invention provides the bivariate nonlocal fuzzy adaptive non-linear average filtering de-noising method for the X-ray image, in the method, a quantum noise model which is hard to process is converted into a common white gaussian noise model, the size of a window of a filter is selected by virtue of fuzzy computation, and a relevant weight matrix enabling an error function to be minimum is searched. A particle swarm optimization filtering parameter is introduced in the method, so that the weight matrix can be locally rebuilt, the influence of the local relevancy on the sample data can be reduced, the algorithm convergence rate can be improved, and the de-noising speed and precision for the industrial X-ray scan image can be improved, so that the method is suitable for processing the X-ray scan image with an uncertain noise model.

Description

The non local average filter radioscopic image of bivariate noise-eliminating method

Technical field

The invention belongs to radioscopic image de-noising field, be specifically related to a kind of noise-eliminating method of radioscopic image of the non local average filter of bivariate (Fuzzy Adaptive Non Local means is abbreviated as FANL means) algorithm of fuzzy self-adaption parameter adjustment.

Background technology

Along with the continuous development of industrial X-ray inspection technique, the quality of X-ray scanning image has also been proposed increasing requirement, this just requires effectively to eliminate the noise information that produces in the real-time testing process.Because narrow, the defective edge fog in X-ray scanning gradation of image interval, many, the covered sometimes characteristics of defect characteristic of picture noise have influenced the effect of workpiece to be detected being analyzed and being estimated according to radioscopic image.Along with the continual renovation of X ray collecting device, new X-ray scanning function is described the information of industrial picture more comprehensively and accurately, and this just causes containing in the image more pixel, and the amount of redundancy between its pixel increases greatly.Therefore, optimize incoherent amount in the operator removal information, and to keep the unchangeability of its inner structure be the gordian technique that non-local average de-noising needs to be resolved hurrily through fuzzy self-adaption.

Non local average filter algorithm utilizes redundant information to carry out denoising Processing; It is an innovation to the local de-noising model of tradition; Its main thought is not that the gray-scale value with single pixel in the image compares; But the distribution situation of the whole gray scale around this pixel is compared, estimate weights according to the similarity of whole intensity profile, its essence is that the relevance map between the original image pixel is handled to image space.It is the de-noising service that Buades points out to make full use of redundant information; In recent years theoretical analysis and the experimental result studied at Buades show; The non local average filter de-noising algorithm of tradition all is superior to common image de-noising algorithm on subjective and objective performance; As gaussian filtering, anisotropy filtering, total error minimize, neighborhood filtering or the like, it originates from the neighborhood filtering algorithm, is a kind of popularization to the neighborhood filtering algorithm; Its weights obtain according to the similarity that whole area grayscale around the pixel distributes, and when reducing picture noise, have the ability of very strong maintenance image spatial resolution.When utilizing traditional non local average filter algorithm process complicated image, its calculated amount is bigger, and processing speed is slow, and especially when handling big image, this problem is more outstanding; In addition, this method can be introduced artificial pseudo-shadow at the smooth region of image, and image thickens, and spatial resolution is affected.

The noise-eliminating method of the radioscopic image of the non local average filter of bivariate of fuzzy self-adaption parameter adjustment is to utilize fuzzy algorithm to realize the self-adaptation adjustment of filtering parameter selection window; Notion and theory with the particle group optimizing operator are the basis; When a plurality of optimization filtering parameters need be selected, utilize fuzzy control rule to carry out the optimum filtering parameters of choice, accomplish evolutionary search through selecting the optimization window to open; Obtain better treatment effect, faster the ability of speed of convergence and global optimizing.But up to the present, also have no talent this algorithm application in non local average optimization aspect.

To the problems referred to above, the present invention improves traditional NL-means method, and algorithm is sought the weighted value that makes reconstruction error minimum from the intrinsic property between data, reaches the purpose that makes error function minimum.Simultaneously, the essence of intelligent optimization algorithm has determined Weight Determination not rely on the distance between data point, reduces the correlativity between the image neighborhood territory pixel effectively; Computational complexity reduces; Saved computing time, improved the algorithm travelling speed, to adapt to the detection needs in the commercial production task.

Summary of the invention

The present invention studies the de-noising situation of the random noise model that produces in the industrial X-ray scan image, adopts the effective denoising Processing that has realized the industrial X-ray scan image based on the quick noise-eliminating method of the non local average filter of the bivariate of fuzzy self-adaption parameter adjustment.

The fuzzy self-adaption parameter adjustment is effectively to the parallel processing of X ray block image; Reduce the correlativity between the image block, realized the quick denoising Processing of X-ray scanning image, simultaneously; Introduce x in the invention; The y bivariate has guaranteed the location invariance of image de-noising process, and, because it adopts particle group optimizing method also can obtain convergence preferably.Thereby the present invention utilizes particle swarm optimization algorithm from the intrinsic property between data, seeks the weighted value that in the X-ray scanning image, makes reconstruction error minimum, thereby makes industrial X-ray de-noising image under the condition that guarantees precision, obtain processing speed faster.

Below the inventive method is made further description, particular content is following:

The non local average filter radioscopic image of bivariate noise-eliminating method, characteristic of the present invention is that method is:

1). the system of selection of fuzzy de-noising window

Non local average filter algorithm has a prerequisite hypothesis: sampled data place local space is linear, and promptly each sampled point and its neighbour point has similarity relation, through the linear expression of weight contribution value;

The learning objective of this algorithm is: in lower dimensional space, make full use of the redundancy relationship between pixel; Similitude based on intensity profile is provided with the weight in each neighborhood; Suppose that promptly the map pane that is bumped into is under the condition of linearity, to minimize uncorrelated pixel, the reconstruct original image in the part;

If c (x y) is the X-ray scanning image function, r (x y) is the ideal image function, n (x y) is the noise image function, and x, y are the coordinate under the rectangular coordinate system of image slices vegetarian refreshments, then have:

c(x，y)＝r(x，y)+n(x，y) (1)

Need look at present a template confirm c (x, y) with r (x, the y) correlativity between each element, closely and effectively eliminate n (x, y);

(the y corresponding gray representes then have with c (i) to c for x, y) pixel x in the image

$\underset{i=1}{\overset{I}{\mathrm{\Σ}}}c\left(i\right)=\underset{y=1}{\overset{Y}{\mathrm{\Σ}}}\underset{x=1}{\overset{X}{\mathrm{\Σ}}}c(x,y)+\mathrm{\γ}---\left(2\right)$

Wherein, I=X*Y is the size dimension of image, and i is picture point x, the gray-scale pixels among the y, i ∈ (1; 2....I), (1,2....X), (1,2....Y) c (i) is the gray-scale value that i is ordered in the X-ray scanning image to y ∈ to x ∈; (half-tone information at y place, γ are the gradation of image modifying factor to c for x, y) reflection image mid point x; Therefore, through regulating every bit γ value, the local contrast of regulating the X-ray scanning image distributes gradation of image after the denoising and suits;

1. random noise model conversion

The main cause of x-ray imaging system image deterioration is system's random noise, on the X-ray scanning image, shows as level or vertical striped, and research shows; The generation of X ray and the noise that produces with the interaction of material; Discrete state and discrete, continuous Markovian process of time are thought in its distribution, all satisfy the Poisson stochastic process in time with on the space, if process with independent increments Δ c (t); The probability distribution of its increment is obeyed Poisson distribution, then has:

$p\left(k\right)\≈\frac{{\mathrm{\λ}}^k}{k!}{e}^{-\mathrm{\λ}}---\left(3\right)$

Wherein, k is the number of times that stochastic variable increment Delta c (t) occurs, k ∈ (1; 2...K); When k is big, according to very inconvenient in the practical application of Poisson distribution formula calculation error ten minutes complicacy, therefore; Adopt among this paper and noise image is carried out this base of a fruit make (Stirling) approximate formula, poisson noise is converted into white Gaussian noise;

Si Di makes (Stirling) approximate formula think, when k is big, Try to achieve by formula (4) is approximate:

$k!\≈\sqrt{2\mathrm{\πk}}{k}^k{e}^{-k}---\left(4\right)$

At this moment, noise is converted into Gaussian distribution from Poisson distribution, and then X ray noise profile expression formula is:

$p\left(k\right)=\frac{1}{\sqrt{2\mathrm{\πh}}}{e}^{[-\frac{{(h-k)}^2}{2h}]}---\left(5\right)$

Then original image meets the white Gaussian noise model, and its mathematical expectation E is asked by formula (6):

$E{\left|\right|c(x_i,y_i)-c(x_j,y_j)\left|\right|}_{2a}^2=E{\left|\right|r(x_i,y_i)-r(x_j,y_j)\left|\right|}_{2a}^2+{2\mathrm{\σ}}^2---\left(6\right)$

Wherein,

is pixel i; J place neighborhood Gauss's weighted euclidean distance; A＞0 is the Gaussian convolution nuclear of standard, and noise variance is σ, the size of the fuzzy de-noising window of size decision of σ;

2. noise level is estimated

Noise level and module window size and filter parameter have confidential relation, therefore, reach denoising effect preferably, just need to estimate the noise level in the piece image.For white Gaussian noise, average is 0, only needs to estimate its standard deviation parameter;

The Gaussian noise method of estimation has a variety of, and commonly used have based on image block with based on two kinds of methods of wave filter.Method based on image block is divided into many fritters with image, calculates its each auto-variance, with the average of some minimum variances as estimated result; Earlier image is carried out once smoothly based on the method for wave filter, calculate the difference of former figure and level and smooth back image again, thus error image estimating noise level.Utilize the piecemeal method that image is divided into many sub-pieces among the present invention earlier, the method for utilizing parallel filtering utilizes particle swarm optimization algorithm to find out maximum noise level in the sub-piece, noise level estimated result as a whole to every sub-block filtering again;

Theoretical analysis that Buades studies in recent years and experimental result point out that the NL-Means algorithm all is superior to common image de-noising algorithm on subjective and objective performance; As gaussian filtering, anisotropy filtering, total error minimize, neighborhood filtering or the like; Its weights obtain according to the similarity that whole area grayscale around the pixel distributes, and when reducing picture noise, have the ability of very strong maintenance image spatial resolution;

3. fuzzy noise filtering template window

In order to improve counting yield; Select two window concurrent operations, through two window sizes are set, one is neighborhood of pixels window size A * A; Another is the window size B * B of neighborhood of pixels window hunting zone; Promptly selecting the neighborhood size of pixel in the inside, zone of B * B size is A * A, carries out self-adaptation non local average (Adaptive Non Local means, ANL-means made in brief note) algorithm; The window of B * B slides in the zone of A * A size, confirms the contribution weights of regional center pixel grey scale according to the similarity in zone.

2). the non local average filter algorithm of bivariate fuzzy self-adaption

Most crucial problem is to utilize the gaussian additive noise level to obtain to make the minimum relevant partial reconstruction weight matrix of reconstruction error in the FANL means algorithm; Yet; This algorithm is to operate to topography's piece, and the researcher adopts the variable relevant with Euclidean distance to define this weight matrix, weighs the size of influence through the distance of distance; This makes this algorithm very sensitive to the noise in the sample, and this algorithm the convergence speed is fast inadequately in addition.The non local average filter algorithm of bivariate fuzzy self-adaption is based upon fuzzy de-noising window and selects on the basis; Each little image block is carried out non local average filter denoising; Make and all obtain a pair of optimization filtering parameter in each module; And utilize particle swarm optimization algorithm to realize the renewal operation of optimized parameter, thereby the optimization that has realized target is found the solution.So the present invention utilizes the optimal parameter of particle swarm optimization algorithm from wave filter, seek and make the minimum weighted value of X-ray scanning characteristics of image vector dimensionality reduction reconstruction error, thereby obtain effective denoising method of X-ray scanning image.

The present invention carries out ANL means Filtering Processing in image subblock; Be configured to the individuality of particle swarm optimization algorithm according to the relative program of sample point and its point of proximity; The individuality of population finds global optimum's speed in the process of optimizing then, confirms the global optimum position, finally obtains the most relevant neighbour's partial reconstruction weight matrix; Thereby remove incoherent redundant information effectively, realize effective de-noising of X ray;

1. the non local average filter of bivariate

If c (x y) is a width of cloth observed image, i.e. X ray detected image, and n is that average is 0, variance is σ ²Additive noise, then (x y) is mapped in the observation space (x, y) two-dimentional bounded domain of definition, (x, y) ∈ R through non local average (Non Local means, NL means made in brief note) algorithm input picture r ², x _i, y _iWith neighborhood point x _j, y _jBetween correlation NL (c (x _i, y _i)) obtain by (6) formula:

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\∫e^{-\frac{(G_a*\left({|c(x_i,y_i)-c(x_j,y_j)|}^2\right)\left(0\right))}{h^2}}c(x_j,y_j)\mathrm{dxdy}---\left(6\right)$

Wherein variable be (x, y),

For standard variance is the Gaussian convolution nuclear of a, h is a filtering parameter, and mainly the noise criteria difference by image determines D (x _i, y _i) be NL means conversion coefficient, it according to coordinate (x, the gray-scale value of relative position y) is obtained:

$D(x_i,y_i)={\∫e}^{\frac{(G_a*\left({\left|c(x_i,y_i)c(x_t,y_t)\right|}^2\right)\left(0\right))}{h^2}}{x}^{\′}\left(i\right)y^{\′}\left(i\right)\mathrm{di}---\left(7\right)$

Be expressed as c (i) for digital X ray scanning image in discrete domain; Its pixel i ∈ (x, y), C={c (i); I ∈ I}; I is a collection of pixels in the image, then the pixel NL means algorithm form that is converted into formula (8) according to the most relevant neighbour's partial reconstruction weight matrix ask between pixel i and the neighborhood correlation NL (c (i)) with relevant neighbour's partial reconstruction weight w (i, j) and the gray-scale value c (j) of its neighbor point represent:

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\∈I}{\mathrm{\Σ}}w(i,j)c\left(j\right)---\left(8\right)$

Wherein, in the formula (2), point (x _i, y _i) gray-scale value located is c (i), (x y) changes correction factor γ and transforms its neighbor point (x by c _j, y _j) gray-scale value be c (j), I is that the pixel of X ray detected image is always counted, X is that (x, width value y), Y are c (x, height value y) to c;

Wherein, i={1,2...i...I}={ (1,1), (1,2) ... (xi, yi) ... (X, Y) }, at this moment, relevant neighbour's partial reconstruction weight w (i j) is expressed as:

$w(i,j)=\frac{1}{Z\left(i\right)}{e}^{{\left|\right|c(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(9\right)$

$Z\left(i\right)=e^{-{\left|\right|c(\mathrm{NL}(x_{i,},y_i)-c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(10\right)$

Wherein, relevant neighbour's partial reconstruction weight w between pixel i and the neighbor point j (i s) sets up the relation of formula (9), and Z (i) is a weight weighting transformation ratio,

Be the Gauss's weighted euclidean distance based on gray level of two neighborhoods at an i and j place, (i s) is mainly determined by parameter a and h weighted value w; A is the standard deviation of the Gaussian convolution nuclear of standard, and h is the filtering parameter of bivariate NL means, and a is big more; Weights are just more little, show pixel x _i, y _iThe distance center pixel is far away more, and the size of a is by the window size A * A decision of selected neighborhood of pixels; Along with the increase of filtering parameter h, artificial pseudo-shadow weakens gradually, but spatial resolution also can descend, and filtering parameter h is by the standard deviation decision of noise;

2. particle group optimizing

Being located at the local optimum parameter that obtains in the p sub-module is a _pAnd h _p, be the local optimum solution vector, utilize the population optimal algorithm, each particle can be estimated the adaptive value of self-position through certain rule; Each particle can be remembered the own current local desired positions z that finds _pPromptly $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right],$ The local location of particle is z _Pq, z _p∈ z _Pq={ z _P1, z _P2...., z _Pq..z _PQ, z _PqBe particle position in the local group kind, Q is the population that comprises in the population, remembers the desired positions that all individualities find in the colony in addition, is called the z of global optimum _g, promptly $\left[\begin{array}{c}a\\ h\end{array}\right],$ Select in the population best $\left[\begin{array}{c}a\\ h\end{array}\right]$ Solution vector;

Particle swarm optimization algorithm points out that all particles all have speed decision their direction and a distance, is called local optimum speed, uses v _pExpression; Particles are followed current optimal particle and in solution space, are searched for, and in each iteration, particle upgrades position and speed according to following formula:

z _p+1＝z _p+v _p+1 (11)

v _p+1＝w ₁v _p+g ₁(s _pq-z _p)+g ₂(s _gq-z _p) (12)

In the formula, Q is the dimension in target search space, and P representes quantity individual in the population, calculates s _pThe current location adaptive value, v _pBe the flying speed of particle p, s _p∈ (s _P1, s _P2, s _Pq... s _PQ) optimal location that searches up to now for particle, s _PqBe the position of particle process in the population, s _pFor speed is v _pThe time optimal location that searches; s _g=(s _G1+ s _G1+ ... s _Gq+ ... s _GQ) optimal location that searches up to now for whole population and, s _GqBe the optimum position of population current record, g ₁And g ₂Be the self study factor, w ₁Be inertia weight;

3. relevance grade evaluation function

The relevance grade evaluation function is to weigh individual good and bad sign, and its effect is in piecemeal subgraph module, to find out whole optimal parameter a and h.The present invention is according to the uncertainty of the industrial X-ray scan image noise model singularity as individual model, and the relevance grade evaluation function that is adopted is:

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2---\left(13\right)$

In the formula, u _PqThe q individuals in p generation in the expression population, D representes quantity individual in the population, the concrete steps of this algorithm optimization are following:

1. initialization population: select threshold epsilon and maximum iteration time u _Max, primary iteration number of times u=0, the particle position scope is z _Min～z _Max, the particle rapidity scope is v _Max～v _Min

2. the individuality in the population is measured: the adaptive value u that measures each particle _Pq, choose this colony's optimal-adaptive value and historical colony optimal-adaptive value compares, obtain the optimal-adaptive value u of colony till current _p, and u _p∈ { u _Pq| u _P1, u _P2... u _Pq... u _PQ, and the value zg ∈ { z of access relevant position _Pq| z _P1, z _P2... z _Pq... z _PQ; Each particle obtains the optimal-adaptive value u till own current through self adaptive value size more in history _g∈ { u _Gq| u _G1, u _G2... u _Gq... u _GQ, and access relevant position z _g∈ { z _Gq| z _G1, z _G2... z _Gq... z _GQ;

3. estimate population U _p, and keep optimum solution;

4. population operation: if U _p＞ε, order is carried out, otherwise end loop;

5. according to rule, upgrade v _p, z _p

6. change evolutionary generation: evolutionary generation adds 1, like the maximum evolutionary generation P of unmet still _MaxAlgorithm goes to this joint step and 2. proceeds.

Beneficial effect of the present invention does; In order better to eliminate the unknown quantum The noise that exists in the industrial X-ray scan image; Proposed to transfer not tractable quantum noise model to common gaussian additive noise model; The size of utilization fuzzy operation selective filter window, searching makes the radioscopic image noise-eliminating method of the non-linear average filter of bivariate fuzzy self-adaption of the minimum relevant weight matrix of error function.In the present invention; Introduce the particle group optimizing filtering parameter; And then the partial reconstruction weight matrix, reduced the influence of local correlations to sample data, improved algorithm the convergence speed; Improved the speed and the precision of the denoising of industrial X-ray scan image, be applicable to the uncertain X-ray scanning treatment of picture of noise model.

Description of drawings

Fig. 1 is the particle swarm optimization algorithm process flow diagram;

Fig. 2 is an industrial X-ray scan image denoising schematic diagram;

Fig. 3 is the subordinate function synoptic diagram of fuzzy window type under the different noise levels.

Embodiment

The non local average filter radioscopic image of bivariate noise-eliminating method, the inventive method do,

1). the system of selection of fuzzy de-noising window

Non local average filter algorithm has a prerequisite hypothesis: sampled data place local space is linear, and promptly each sampled point and its neighbour point has similarity relation, through the linear expression of weight contribution value;

The learning objective of this algorithm is: in lower dimensional space, make full use of the redundancy relationship between pixel; Similitude based on intensity profile is provided with the weight in each neighborhood; Suppose that promptly the map pane that is bumped into is under the condition of linearity, to minimize uncorrelated pixel, the reconstruct original image in the part;

If c (x y) is the X-ray scanning image function, r (x y) is the ideal image function, n (x y) is the noise image function, and x, y are the coordinate under the rectangular coordinate system of image slices vegetarian refreshments, then have:

c(x，y)＝r(x，y)+n(x，y) (1)

Need look at present a template confirm c (x, y) with r (x, the y) correlativity between each element, closely and effectively eliminate n (x, y);

(the y corresponding gray representes then have with c (i) to c for x, y) pixel x in the image

$\underset{i=1}{\overset{I}{\mathrm{\Σ}}}c\left(i\right)=\underset{y=1}{\overset{Y}{\mathrm{\Σ}}}\underset{x=1}{\overset{X}{\mathrm{\Σ}}}c(x,y)+\mathrm{\γ}---\left(2\right)$

Wherein, I=X*Y is the size dimension of image, and i is picture point x, the gray-scale pixels among the y, i ∈ (1; 2....I), (1,2....X), (1,2.....Y) c (i) is the gray-scale value that i is ordered in the X-ray scanning image to y ∈ to x ∈; (half-tone information at y place, γ are the gradation of image modifying factor to c for x, y) reflection image mid point x; Therefore, through regulating every bit γ value, the local contrast of regulating the X-ray scanning image distributes gradation of image after the denoising and suits;

1. random noise model conversion

The main cause of x-ray imaging system image deterioration is system's random noise, on the X-ray scanning image, shows as level or vertical striped, and research shows; The generation of X ray and the noise that produces with the interaction of material; Discrete state and discrete, continuous Markovian process of time are thought in its distribution, all satisfy the Poisson stochastic process in time with on the space, if process with independent increments Δ c (t); The probability distribution of its increment is obeyed Poisson distribution, then has:

$p\left(k\right)\≈\frac{{\mathrm{\λ}}^k}{k!}{e}^{-\mathrm{\λ}}---\left(3\right)$

Wherein, k is the number of times that stochastic variable increment Delta c (t) occurs, k ∈ (1; 2...K); When k is big, according to very inconvenient in the practical application of Poisson distribution formula calculation error ten minutes complicacy, therefore; Adopt among this paper and noise image is carried out this base of a fruit make (Stirling) approximate formula, poisson noise is converted into white Gaussian noise;

Si Di makes (Stirling) approximate formula think, when k is big, Try to achieve by formula (4) is approximate:

$k!\≈\sqrt{2\mathrm{\πk}}{k}^k{e}^{-k}---\left(4\right)$

At this moment, noise is converted into Gaussian distribution from Poisson distribution, and then X ray noise profile expression formula is:

$p\left(k\right)=\frac{1}{\sqrt{2\mathrm{\πh}}}{e}^{[-\frac{{(h-k)}^2}{2h}]}---\left(5\right)$

Then original image meets the white Gaussian noise model, and its mathematical expectation E is asked by formula (6):

$E{\left|\right|c(x_i,y_i)-c(x_j,y_j)\left|\right|}_{2a}^2=E{\left|\right|r(x_i,y_i)-r(x_j,y_j)\left|\right|}_{2a}^2+{2\mathrm{\σ}}^2---\left(6\right)$

Wherein,

is pixel i; J place neighborhood Gauss's weighted euclidean distance; A＞0 is the Gaussian convolution nuclear of standard, and noise variance is σ, the size of the fuzzy de-noising window of size decision of σ;

2. noise level is estimated

Noise level and module window size and filter parameter have confidential relation, therefore, reach denoising effect preferably, just need to estimate the noise level in the piece image.For white Gaussian noise, average is 0, only needs to estimate its standard deviation parameter;

The Gaussian noise method of estimation has a variety of, and commonly used have based on image block with based on two kinds of methods of wave filter.Method based on image block is divided into many fritters with image, calculates its each auto-variance, with the average of some minimum variances as estimated result; Earlier image is carried out once smoothly based on the method for wave filter, calculate the difference of former figure and level and smooth back image again, thus error image estimating noise level.Utilize the piecemeal method that image is divided into many sub-pieces among the present invention earlier, the method for utilizing parallel filtering utilizes particle swarm optimization algorithm to find out maximum noise level in the sub-piece, noise level estimated result as a whole to every sub-block filtering again;

Theoretical analysis that Buades studies in recent years and experimental result point out that the NL-Means algorithm all is superior to common image de-noising algorithm on subjective and objective performance; As gaussian filtering, anisotropy filtering, total error minimize, neighborhood filtering or the like; Its weights obtain according to the similarity that whole area grayscale around the pixel distributes, and when reducing picture noise, have the ability of very strong maintenance image spatial resolution;

3. fuzzy noise filtering template window

In order to improve counting yield; Select two window concurrent operations, through two window sizes are set, one is neighborhood of pixels window size A * A; Another is the window size B * B of neighborhood of pixels window hunting zone; Promptly selecting the neighborhood size of pixel in the inside, zone of B * B size is A * A, carries out self-adaptation non local average (Adaptive Non Local means, ANL-means made in brief note) algorithm; The window of B * B slides in the zone of A * A size, confirms the contribution weights of regional center pixel grey scale according to the similarity in zone.

2). the non local average filter algorithm of bivariate fuzzy self-adaption

Most crucial problem is to utilize the gaussian additive noise level to obtain to make the minimum relevant partial reconstruction weight matrix of reconstruction error in the FANL means algorithm; Yet; This algorithm is to operate to topography's piece, and the researcher adopts the variable relevant with Euclidean distance to define this weight matrix, weighs the size of influence through the distance of distance; This makes this algorithm very sensitive to the noise in the sample, and this algorithm the convergence speed is fast inadequately in addition.The non local average filter algorithm of bivariate fuzzy self-adaption is based upon fuzzy de-noising window and selects on the basis; Each little image block is carried out non local average filter denoising; Make and all obtain a pair of optimization filtering parameter in each module; And utilize particle swarm optimization algorithm to realize the renewal operation of optimized parameter, thereby the optimization that has realized target is found the solution.So the present invention utilizes the optimal parameter of particle swarm optimization algorithm from wave filter, seek and make the minimum weighted value of X-ray scanning characteristics of image vector dimensionality reduction reconstruction error, thereby obtain effective denoising method of X-ray scanning image.

The present invention carries out ANL means Filtering Processing in image subblock; Be configured to the individuality of particle swarm optimization algorithm according to the relative program of sample point and its point of proximity; The individuality of population finds global optimum's speed in the process of optimizing then, confirms the global optimum position, finally obtains the most relevant neighbour's partial reconstruction weight matrix; Thereby remove incoherent redundant information effectively, realize effective de-noising of X ray;

1. the non local average filter of bivariate

If c (x y) is a width of cloth observed image, i.e. X ray detected image, and n is that average is 0, variance is σ ²Additive noise, then (x y) is mapped in the observation space (x, y) two-dimentional bounded domain of definition, (x, y) ∈ R through non local average (Non Local means, NL means made in brief note) algorithm input picture r ², x _i, y _iWith neighborhood point x _j, y _jBetween correlation NL (c (x _i, y _i)) obtain by (6) formula:

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\∫e^{-\frac{(G_a*\left({|c(x_i,y_i)-c(x_j,y_j)|}^2\right)\left(0\right))}{h^2}}c(x_j,y_j)\mathrm{dxdy}---\left(6\right)$

Wherein variable be (x, y),

For standard variance is the Gaussian convolution nuclear of a, h is a filtering parameter, and mainly the noise criteria difference by image determines D (x _i, y _i) be NL means conversion coefficient, it according to coordinate (x, the gray-scale value of relative position y) is obtained:

$D(x_i,y_i)={\∫e}^{\frac{(G_a*\left({\left|c(x_i,y_i)c(x_t,y_t)\right|}^2\right)\left(0\right))}{h^2}}{x}^{\′}\left(i\right)y^{\′}\left(i\right)\mathrm{di}---\left(7\right)$

Be expressed as c (i) for digital X ray scanning image in discrete domain; Its pixel i ∈ (x, y), C={c (i); I ∈ I}; I is a collection of pixels in the image, then the pixel NL means algorithm form that is converted into formula (8) according to the most relevant neighbour's partial reconstruction weight matrix ask between pixel i and the neighborhood correlation NL (c (i)) with relevant neighbour's partial reconstruction weight w (i, j) and the gray-scale value c (j) of its neighbor point represent:

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\∈I}{\mathrm{\Σ}}w(i,j)c\left(j\right)---\left(8\right)$

Wherein, in the formula (2), point (x _i, y _i) gray-scale value located is c (i), (x y) changes correction factor γ and transforms its neighbor point (x by c _j, y _j) gray-scale value be c (j), I is that the pixel of X ray detected image is always counted, X is that (x, width value y), Y are c (x, height value y) to c;

Wherein, i={1,2...i...I}={ (1,1), (1,2) ... (x _i, y _i) .... (X, Y) }, at this moment, relevant neighbour's partial reconstruction weight w (i j) is expressed as:

$w(i,j)=\frac{1}{Z\left(i\right)}{e}^{{\left|\right|c(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(9\right)$

$Z\left(i\right)=e^{-{\left|\right|c(\mathrm{NL}(x_{i,},y_i)-c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(10\right)$

Wherein, relevant neighbour's partial reconstruction weight w between pixel i and the neighbor point j (i s) sets up the relation of formula (9), and Z (i) is a weight weighting transformation ratio, Be the Gauss's weighted euclidean distance based on gray level of two neighborhoods at an i and j place, (i s) is mainly determined by parameter a and h weighted value w; A is the standard deviation of the Gaussian convolution nuclear of standard, and h is the filtering parameter of bivariate NL means, and a is big more; Weights are just more little, show pixel x _i, y _iThe distance center pixel is far away more, and the size of a is by the window size A * A decision of selected neighborhood of pixels; Along with the increase of filtering parameter h, artificial pseudo-shadow weakens gradually, but spatial resolution also can descend, and filtering parameter h is by the standard deviation decision of noise;

2. particle group optimizing

Being located at the local optimum parameter that obtains in the p sub-module is a _pAnd h _p, be the local optimum solution vector, utilize the population optimal algorithm, each particle can be estimated the adaptive value of self-position through certain rule; Each particle can be remembered the own current local desired positions z that finds _pPromptly $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right],$ The local location of particle is z _Pq, z _p∈ z _Pq={ z _P1, z _P2...., z _Pq..z _PQ, z _PqBe particle position in the local group kind, Q is the population that comprises in the population, remembers the desired positions that all individualities find in the colony in addition, is called the z of global optimum _g, promptly $\left[\begin{array}{c}a\\ h\end{array}\right],$ Select in the population best $\left[\begin{array}{c}a\\ h\end{array}\right]$ Solution vector;

Particle swarm optimization algorithm points out that all particles all have speed decision their direction and a distance, is called local optimum speed, uses v _pExpression; Particles are followed current optimal particle and in solution space, are searched for, and in each iteration, particle upgrades position and speed according to following formula:

z _p+1＝z _p+v _p+1 (11)

v _p+1＝w ₁v _p+g ₁(s _pq-z _p)+g ₂(s _gq-z _p) (12)

In the formula, Q is the dimension in target search space, and P representes quantity individual in the population, calculates s _pThe current location adaptive value, v _pBe the flying speed of particle p, s _p∈ (s _P1, s _P2, s _Pq... s _PQ) optimal location that searches up to now for particle, s _PqBe the position of particle process in the population, s _pFor speed is v _pThe time optimal location that searches; s _g=(s _G1+ s _G1+ ... s _Gq+ ... s _GQ) optimal location that searches up to now for whole population and, s _GqBe the optimum position of population current record, g ₁And g ₂Be the self study factor, w ₁Be inertia weight;

3. relevance grade evaluation function

The relevance grade evaluation function is to weigh individual good and bad sign, and its effect is in piecemeal subgraph module, to find out whole optimal parameter a and h.The present invention is according to the uncertainty of the industrial X-ray scan image noise model singularity as individual model, and the relevance grade evaluation function that is adopted is:

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2---\left(13\right)$

In the formula, u _PqThe q individuals in p generation in the expression population, D representes quantity individual in the population, the concrete steps of this algorithm optimization are following:

1. initialization population: select threshold epsilon and maximum iteration time u _Max, primary iteration number of times u=0, the particle position scope is z _Min～z _Max, the particle rapidity scope is v _Max～v _Min

2. the individuality in the population is measured: the adaptive value u that measures each particle _Pq, choose this colony's optimal-adaptive value and historical colony optimal-adaptive value compares, obtain the optimal-adaptive value u of colony till current _p, and u _p∈ { u _Pq| u _P1, u _P2... u _Pq... u _PQ, and the value z of access relevant position _g∈ { z _Pq| z _P1, z _P2... z _Pq... z _PQ; Each particle obtains the optimal-adaptive value u till own current through self adaptive value size more in history _g∈ { u _Gq| u _G1, u _G2... u _Gq... u _GQ, and access relevant position z _g∈ { z _Gq| z _G1, z _G2... z _Gq... z _GQ;

3. estimate population U _p, and keep optimum solution;

4. population operation: if U _p＞ε, order is carried out, otherwise end loop;

5. according to rule, upgrade v _p, z _p

6. change evolutionary generation: evolutionary generation adds 1, like the maximum evolutionary generation P of unmet still _MaxAlgorithm goes to this joint step and 2. proceeds.

The characteristics such as X-ray scanning gradation of image interval is narrow, defective edge fog, picture noise is many, defect characteristic is submerged sometimes that the present invention is directed in the electric field commercial unit Non-Destructive Testing process are carried out Denoising Study; Respectively with horizontal quantum noise; Vertical quantum noise; Random noise, the noise of crossing the enhancement process generation is that example is studied, multiple format-patterns such as input BMP, JPG, PNG are the scan image of size arbitrarily.At first according to the denoising index request auto-selecting parameter (target image) is proposed by the user; Do not provide then like the user and to differentiate automatically, and convert the quantum noise of the image of input into Gaussian noise, then by system; High phase noise is carried out grade classification; Distribution according to the noise model level utilizes fuzzy control rule that different windows is set, and in each subwindow, the different filtering parameter is set, and utilizes the parallel particle colony optimization algorithm to carry out the selection of global optimum's filtering parameter; Closely view picture X-ray scanning image is carried out the FANL Denoising Algorithm; And with least mean-square error (Minimum Mean SquaredError, i.e. MMSE), Y-PSNR (peak signal noise ratio, i.e. PSNR) and the used time t of denoising evaluation system denoising performance.System architecture diagram is as shown in Figure 2, and the practical implementation method is following:

1. noise level modeling and template window size Selection

The main cause of x-ray imaging system image deterioration is system's random noise.Research shows, the generation of X ray and with the interaction of material, all satisfy the Poisson stochastic process in time with on the space.For the rapid X-ray imaging system, because the time shutter is short, the quantum noise that X ray produced is more outstanding, has had a strong impact on the quality of image.Therefore, to the X-ray scanning image of input,, therefore, on the original image basis, add average μ because its noise model is unknown ₀Be 0, initial variance σ ₀Be 0.02 low noise level at white Gaussian noise, wherein, μ ₀, σ ₀For adding the initial value of Gaussian noise; It is a normal value; At this moment, utilize this base of a fruit to make (Stirling) formula (4) that the image function model that X-ray scanning adds behind the even noise of Gauss is similar to, poisson noise is converted into white Gaussian noise; Obtain the variances sigma of the overall white Gaussian noise of conversion back image again, the σ value influences system performance.

Select principle for use according to fuzzy noise filtering template window, the image that need add the radioscopic image after detecting after Gauss's noise processed blur each local optimum parameter of block parallel calculating, therefore; Two windows need be set, neighborhood of pixels window size A * A, the window B * B of neighborhood window hunting zone; The image that common noise is bigger is generally got A=7, B=23; Get A=3 for low noise image, B=7 just can satisfy noise reduction basically.Therefore, noise level is provided with the template window size of non local average de-noising, and adding parameter earlier for the noise of importing is (u ₀, σ ₀) white Gaussian noise after the conversion noise model; Difference according to the high phase noise level σ in the image function of conversion back is divided into five grades; Be low noise σ＜0.2, than low noise 0.2＜σ＜0.5, medium noise 0.5＜σ＜0.8; Higher noise 0.8＜σ＜1 and a domain interval, strong noise σ＞1 five; For different noises different processing template windows is set, in each template window, adjusts filtering parameter adaptively, the type level of obscures window as shown in Figure 3 under the varying level.

2. the local optimum parameter is selected in the template window

(1) local parameter a

The weighted value w of pixel i and neighborhood pixels point j in the non local average filter (i, j) mainly by parameter a and h decision, a is the Gaussian convolution nuclear of standard, and h is the filtering parameter of bivariate NL means, and a is big more for Gaussian convolution nuclear, and weights are just more little, show pixel x _i, y _iThe distance center pixel is far away more, and the size of a is by the window size A * A decision of selected neighborhood of pixels; The level and smooth degree of the value of a and image itself also has certain relation; When a value is big; Comprised more noise spot; But also can point that should be different with the current pixel value in the image be comprised to come in simultaneously, so a value is not to be the bigger the better, a often gets in 1～10 times the scope of neighborhood window size; Take into account the influence of different noise levels simultaneously, get 10～15 times of noise criteria difference σ values as a reference, consider that all and current pixel weighting absolute difference are no more than the pixel of a σ.Therefore, when the optimized parameter of selecting is discontented when being enough to condition, then give up the local optimum parameter value of having selected.

(2) local parameter h

Research draws, filtering parameter h and noise variance σ ²The approximately linear proportional relation is arranged, and receives the influence of noise image variance, in order not lose the information outside the noise, when the Weighted distance of two pixels greater than threshold value h _βThe time (β is the cut-off condition of filter configuration), (i j) should approach zero to weight w.

When the distribution of pixel meets Gaussian distribution in adding the processing back image of making an uproar, and its noise weights function w (i, j) also meeting average is 0; Standard deviation is the typical Gaussian function of h/2, therefore, regulates weight w (i through changing the h value; J) distribution, then (i j) satisfies the characteristics of Gaussian distribution to w; The point identical with current pixel point and neighborhood territory pixel gray scale thereof can have some differences with current pixel value after the Gaussian noise (σ) that has superposeed; These differences have determined i, and the j point-to-point transmission is based on Gauss's weighted euclidean distance value of gray level, and the weighting absolute difference that has a certain proportion of pixel to drop on current considered pixel value is no more than within the noise criteria difference β scope doubly.Setting is ignored weights and is ignored noise pixel number percent when calculating Weighted distance, when weights are estimated out apart from all weights sums of a certain multiple of overgauge deviation, set the weights ignored with, ask to obtain this multiple value.The multiple that the weighting absolute difference of at present supposing pixel value is no more than the noise criteria difference is β, and then not uncared-for distance should satisfy:

${\left|\right|c\left({\mathrm{NL}}_i\right)-c\left({\mathrm{NL}}_j\right)\left|\right|}_{2a}^2\≤\frac{\mathrm{\βh}}{\sqrt{2}}---\left(14\right)$

The approximation relation that then gets filtering parameter and noise level is:

$h=\sqrt{2}{\mathrm{\β\σ}}^2---\left(15\right)$

The radioscopic image number of greyscale levels is 256, therefore, adopts following formula to estimate the h value:

$h=[\sqrt{2}\mathrm{\β}+\frac{\max \left(O({\mathrm{\σ}}_c-{\mathrm{\σ}}_0)\right)}{10}]{\mathrm{\σ}}^2---\left(16\right)$

σ _cFor the average of detected image poor, σ ₀Be normative reference difference limen value, when the parameter in the employing formula (16) is carried out filtering to noise image, reach the effect of near-optimization.The too small meeting of h parameter value is not considered many noises; The h parameter value is excessive, can smoothly fall the image difference that some exceed the noise scope, thus cause image too level and smooth, lose detailed information such as edge.For gray scale is 256 grades X-ray scanning image, less than σ ₀Standard deviation less to h value influence, and be higher than σ when noise level ₀In the time, influence bigger to h.

3. select based on the filtering parameter solution vector of particle group optimizing

When radioscopic image is carried out denoising Processing, at first to carry out piecemeal, utilize ANL means algorithm to obtain the local optimum parametric solution vector of every sub-block again entire image.Because the noise model during X-ray production apparatus scanning obtains is unknown; Therefore different to the number of different x-ray image block; Local optimum parameter in each image subblock also is not quite similar, and therefore, adopts particle swarm optimization algorithm in the X-ray scanning image, to seek global optimum's parametric solution vector.

Come characteristic information extraction according to the correlativity of each pixel gray-scale value in the X-ray scanning image subblock of the present invention after the Gaussian noise model conversion, wherein, neighborhood window size A * A; Window size B * the B of hunting zone at this moment, not only is correlated with between the consecutive point; Non-conterminously in setting range, also exist certain correlativity, (i is j) according to pixel i for weighted value w; Correlativity between the j obtains, and in every sub-block, shows the local optimum parameter a that obtains every sub-block _i, h _i, utilize particle swarm optimization algorithm to obtain overall optimal parameter a, h again.

If the input radioscopic image is divided into the P sub-block, consider the correlativity between pixel, also exist certain correlativity between every sub-block, therefore, each sub-block is regarded a sample point c in optimized parameter to be selected space as _p, calculate each sample point (proper vector) c _pK nearest neighbor point, calculate the distance between this point and other P-1 the sample point, will be apart from ordering, k and c before selecting _pNearest point is as its neighbor point.

Calculate the partial reconstruction weight matrix of this sample point by neighbour's point of each sample point $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right].$ Neighbour's point with each proper vector is rebuild this proper vector, asks for the neighbour's partial reconstruction weight matrix that makes reconstruction error minimum.

It is following that the error of fitting function is rebuild in definition:

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2$

Wherein, c _q(q=1,2 ..., k) be c _pQ in the sample neighbour's point, u _PqBe c _pWith c _qBetween weights.

When satisfying following two constraint conditions, obtain the partial reconstruction weight matrix through the minimize error function, promptly by neighbour's point of sample point, construct optimum W matrix and make the error function value reach minimum.

1) each data point c _qAll can only represent, if c by its reference point _qNot reference point, then u _Pq=0;

2) each row of weight matrix and be 1, promptly satisfy normalization constraint

In the process of asking for optimum U matrix; U is for rebuilding the error function that error of fitting function solution vector constitutes; Reconstruction weight vector with sample point and its reference point among the present invention is gathered the primary as particle cluster algorithm; Individuality in the population finds global optimum position and optimum velocity in the process of optimizing then, finally obtains neighbour's partial reconstruction weight matrix U.Wherein, u _PqThe q individuals in p generation in the expression population, P representes quantity individual in the population.The concrete steps of this algorithm optimization are as follows, and its algorithm flow is as shown in Figure 1.

In the formula, u _PqThe q individuals in p generation in the expression population, D representes quantity individual in the population, and the concrete steps of this algorithm optimization are as follows, and its algorithm flow is as shown in Figure 1.

1) initialization population: select threshold epsilon and maximum iteration time u _Max, primary iteration number of times u=0, the particle position scope is z _Min～z _Max, the particle rapidity scope is v _Max～v _Min

2) individuality in the population is measured: the adaptive value u that measures each particle _Pq, choose this colony's optimal-adaptive value and historical colony optimal-adaptive value compares, obtain the optimal-adaptive value u of colony till current _p, and u _p∈ { u _Pq| u _P1, u _P2... u _Pq... u _PQ, and the value z of access relevant position _g∈ { z _Pq| z _P1, z _P2... z _Pq... z _PQ; Each particle obtains the optimal-adaptive value u till own current through self adaptive value size more in history _g∈ { u _Gq| u _G1, u _G2... u _Gq... u _GQ, and access relevant position z _g∈ { z _Gq| z _G1, z _G2... z _Gq... z _GQ.

3) estimate population U _p, and keep optimum solution;

4) population operation: if U _p＞ε, order is carried out, otherwise end loop;

5), upgrade v according to rule _p, z _p

6) change evolutionary generation: evolutionary generation adds 1, like the maximum evolutionary generation P of unmet still _MaxAlgorithm goes to step 2) proceed.

Promptly ask for optimum neighbour's partial reconstruction weight matrix U according to above step, for keeping the weights U (P) that asks for of step constant, algorithm flow the 2nd) constraint condition carried out of step is:

1) location update operations: z _P+1=z _p+ v _P+1

2) Velocity Updating operation: v _P+1=w ₁v _p+ g ₁(s _Pq-z _p)+g ₂(s _Gq-z _p)

In the formula, q ∈ { q ₁, q ₂... Q}, Q are the dimensions in target search space, and P is the number of data point in the sample, calculate s _pThe current location adaptive value, v _pBe the flying speed of particle p, s _p∈ (s _P1, s _P2, s _Pq... s _PQ) optimal location that searches up to now for particle, s _PqBe the position of particle process in the population, s _pFor speed is v _pThe time optimal location that searches; s _g=(s _G1+ s _G1+ ... s _Gq+ ... s _GQ) optimal location that searches up to now for whole population and, s _GqBe the optimum position of population current record, g ₁And g ₂Be the self study factor, w ₁Be inertia weight.Therefore, the corresponding vector of Q eigenwert that minimizes the minimum of getting f (u) in the error of fitting function is that column vector is formed matrix [a _p, h _p] ^T, then the column vector of U is the vector representation of Q dimension space.

4. similarity measurement

Convert the quantum noise in the industrial X-ray scan image of input into Gaussian noise among the present invention, closely the imagery exploitation FANL algorithm that has Gaussian noise after the conversion is carried out quick denoising.Therefore, the radioscopic image that has Gaussian noise after industrial X-ray detected image that has uncertain quantum noise that detects and the noise model conversion that adds low noise level is carried out similarity measurement, the present invention adopts the following similar method of estimating:

(1) mutual information tolerance

Mutual information when the value of mutual information is big more, representes that two magnitude image are similar more as the similarity measure function of weighing two width of cloth image registrations.

Suppose that the industrial X-ray detected image that has uncertain quantum noise that detects is image C R, the radioscopic image that has Gaussian noise after the conversion is an image C, and C is C by selecting the sub collective drawing behind the neighborhood window piecemeal ₁, C ₂..C _a..C _A, C is divided into C by the search neighborhood window ₁, C ₂..C _b.C _B, CR, C are any real number, CR then, and the marginal probability density of C is respectively P _CR(cr) and P _C(c), their joint probability distribution then is expressed as P _AB(a b), then uses the degree of correlation between two signals of simple crosscorrelation information representation, then has

$I(\mathrm{CR},C)=\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log \frac{P_{\mathrm{CRC}}(\mathrm{cr},c)}{P_{\mathrm{CR}}\left(\mathrm{cr}\right)P_C\left(c\right)}---\left(17\right)$

Simple crosscorrelation and entropy have confidential relation, and formula (17) is write the form of an accepted way of doing sth (18):

I(CR，C)＝H(CR，C)＝H(CR)+H(C)-H(CR，C)＝H(C)H(C/CR) (18)

H (A) and H (B) are respectively the entropy of two width of cloth images.Relation table before the information entropy of piece image CR, it and each pixel be shown as shown in the formula form:

$H\left(\mathrm{CR}\right)=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{pi}\log \mathrm{pi}---\left(19\right)$

I ∈ (1,2 ... the .n) pixel count that is contained of expression signal A, establish H (CR, C) presentation video CR, the combination entropy of C, H (CR/C), H (C/CR) be the conditional entropy of presentation video respectively, then has:

$H\left(\mathrm{CR}\right|C)=-\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log P_{\mathrm{CRC}}\left(\mathrm{cr}\right|c)---\left(20\right)$

$H(\mathrm{CR},C)=-\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log P_{\mathrm{CRC}}(\mathrm{cr},c)---\left(21\right)$

If m (e) be the gray-scale value of former detected image at a gray-scale value at s place, e is the pixel in the image, Γ (T (e)) is the gray-scale value of original image dry acoustic model conversion back after an e place carries out the Γ conversion; The gray-scale value size is set between 0～255, and in like manner, Γ (e) is 0～255 with the gray-scale value scope of m (T (e)); When C moves on CR; Lap forms two-dimensional histogram, and the gray scale that note work two width of cloth figure form is h to the number of (Γ (e), m (T (e))) _α(Γ m), then the image C R that is associated of two width of cloth and the joint distribution probability of C do

The marginal distribution probability

With

So cross correlation function is:

$I(\∂)=\underset{\mathrm{cr},c}{\mathrm{\Σ}}{P}_{\mathrm{CRC},\∂}(\mathrm{cr},c){\log }_2\frac{P_{\mathrm{CRC},\∂}(\mathrm{cr},c)}{P_{\mathrm{CR},\∂}\left(\mathrm{cr}\right)P_C\left(c\right)}---\left(22\right)$

Optimum mutual information parameter

has

${\∂}^{*}=\mathrm{arc}\max I\left(\mathrm{\α}\right)---\left(23\right)$

Obtain neighborhood window C respectively ₁, C ₂..C _a..C _AWith search window C ₁, C ₂..C _b.C _BOptimum mutual information parameter

With

Optimal value then ${{\∂}_o}^{*}=\mathrm{Max}\{{\∂}_a^{*},{\∂}_b^{*}\}.$

(2) cosine distance metric

In order to investigate the degree of correlation between each location of pixels of two magnitude image, utilize the method for histogrammic friendship complementation chordal distance method to find out the correlationship on the picture position.

Suppose that F and Q are the histograms that comprises w point of image C and CR image, then crossing between them represented with following formula apart from l:

$l=\underset{w=1}{\overset{N}{\mathrm{\Σ}}}\min (F_w,Q_w)---\left(24\right)$

Histogrammic intersecting is meant two histograms total pixel quantity in each point.Sometimes, this value is also through realizing standardization divided by pixel quantities all in one of them histogram, thereby makes its value be in the codomain scope of [0,1], then has:

$l(F,Q)=\underset{w=1}{\overset{N}{\mathrm{\Σ}}}\min (F_w,Q_w)/\underset{w=1}{\overset{N}{\mathrm{\Σ}}}{Q}_w---\left(25\right)$

Trying to achieve the cosine distance is:

l(F，Q)＝F ^T*Q/(||F||*||Q||) (26)

Wherein, F and Q represent the proper vector of image in query image and the database respectively, || * | the expression vector norm.The similarity measurement value that calculates is between [0,1], and this value is big more, and presentation video is similar more.

For the assemblage characteristic image, similarity measurement is defined as the weighted sum of each characteristic similarity tolerance.Its formula is:

$l(F,Q)=\underset{w=1}{\overset{m}{\mathrm{\Σ}}}{u}_w*l_w(F,Q)---\left(27\right)$

Expression is formed by m characteristics combination, wherein u _wThe weight coefficient of representing w characteristic, it representes the importance of w characteristic, generally gets each u _wEquate S _w(F, Q) the similarity measurement functional value of w characteristic of expression.

5. the performance evaluation of algorithm

The present invention adopts least mean-square error MMSE (Minimum Means SquaredError) commonly used in the denoising performance evaluation, and Y-PSNR PSNR and response time t come the quality of objective evaluation denoising performance.

The intensity of variation of view data after the MMSE evaluation denoising, the value of MMSE is more little, explains that denoising effect is accurate more.The noise model of former X ray detected image is a quantum noise, is σ with the noise square error in the Gaussian noise model that draws after its conversion ₀, image through the method denoising Processing among the present invention after, can know the information entropy of new images each point pixel after the de-noising, with the Pixel Information entropy of the pixel entropy of image after the de-noising and former X-ray scanning image each point relatively, obtain the meansquaredeviation of integral body ₁, MMSE=min{ σ then ₀, σ ₁, the situation of the average de-noising level of image each point after the denoising in the MMSE reflection.For a little bigger influence of difference in the key diagram picture to denoising effect; Introduce PSNR among the present invention as another objective standard of estimating the image denoising effect quality; Be used for weighing the situation of change between image maximum of points and the noise; It is CR that known detection obtains the X-ray scanning image, and it is C that CR is converted to image through noise model, obtains image R after the C de-noising ^*, ideal situation expects that down the not noisy that obtains is R at image, then the PSNR definition is as follows:

$\mathrm{PSNR}=10\×\log \left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(28\right)$

$\mathrm{MMSE}=\min \left\{\right|\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}({\left(R^x\right)}^n-R^n)}{N}\left|\right\}---\left(29\right)$

Wherein, N is a picture size, and n is the pixel in the image, contained pixel count in the n ∈ N presentation video.PSNR and MMSE show the performance index of the effect of this denoising system, according to user's request, the PSNR and the MMSE threshold value of acceptance are set in the reality.

Response time also is an index of ignoring in weighing the denoising system performance index; Among the present invention owing to adopt each module parallel computation behind the piecemeal; Thereby be superior to traditional NL means Denoising Algorithm in the processing time, but consider precision requirement aspect, utilize x (n) among the present invention; Y (n) describes the positional information of pixel n in image, thereby can handle the denoising situation of appointed area more accurately.

When using traditional non local average filter de-noising algorithm process complicated image, its calculated amount is bigger, and processing speed is slow, and especially when handling big image, this problem is more outstanding; In addition, traditional non local average filter method can be introduced artificial pseudo-shadow at the smooth region of image, and image thickens, and spatial resolution is affected.The non local average filter de-noising of the fuzzy bivariate algorithm that proposes among the present invention utilizes the mode of piecemeal to carry out piecemeal to complicated image; The antithetical phrase piece carries out the auto-adaptive parameter concurrent operation again; When reducing image complexity, saved computing time, reached denoising performance preferably.

Claims (1)

1. the non local average filter radioscopic image of bivariate noise-eliminating method is characterized in that, method is:

1). the system of selection of fuzzy de-noising window

Non local average filter algorithm has a prerequisite hypothesis: sampled data place local space is linear, and promptly each sampled point and its neighbour point has similarity relation, through the linear expression of weight contribution value;

This algorithm makes full use of the redundancy relationship between pixel in lower dimensional space, according to the similarity of intensity profile the weight in each neighborhood is set, and promptly the map pane that is bumped into of hypothesis minimizes uncorrelated pixel, the reconstruct original image under the part is linear condition;

If c (x y) is the X-ray scanning image function, r (x y) is the ideal image function, n (x y) is the noise image function, and x, y are the coordinate under the rectangular coordinate system of image slices vegetarian refreshments, then have:

c(x，y)＝r(x，y)+n(x，y) (1)

Need look at present a template confirm c (x, y) with r (x, the y) correlativity between each element, closely and effectively eliminate n (x, y);

(the y corresponding gray representes then have with c (i) to c for x, y) pixel x in the image

$\underset{i=1}{\overset{I}{\mathrm{\Σ}}}c\left(i\right)=\underset{y=1}{\overset{Y}{\mathrm{\Σ}}}\underset{x=1}{\overset{X}{\mathrm{\Σ}}}c(x,y)+\mathrm{\γ}---\left(2\right)$

Wherein, I=X*Y is the size dimension of image, and i is picture point x, the gray-scale pixels among the y, i ∈ (1; 2...I), (1,2...X), (1,2.....Y) c (i) is the gray-scale value that i is ordered in the X-ray scanning image to y ∈ to x ∈; (half-tone information at y place, γ are the gradation of image modifying factor to c for x, y) reflection image mid point x; Therefore, through regulating every bit γ value, the local contrast of regulating the X-ray scanning image distributes gradation of image after the denoising and suits;

1. random noise model conversion

The main cause of x-ray imaging system image deterioration is system's random noise; On the X-ray scanning image, show as level or vertical striped, the generation of X ray and the noise that produces with the interaction of material, its distribution is a continuous Markovian process of discrete state and discrete, time; All satisfy the Poisson stochastic process in time with on the space; If process with independent increments Δ c (t), the probability distribution of its increment is obeyed Poisson distribution, then has:

$p\left(k\right)\≈\frac{{\mathrm{\λ}}^k}{k!}{e}^{-\mathrm{\λ}}---\left(3\right)$

Wherein, k is the number of times that stochastic variable increment Delta c (t) occurs, k ∈ (1; 2...K); When k is big, according to very inconvenient in the practical application of Poisson distribution formula calculation error ten minutes complicacy, therefore; Employing is carried out this base of a fruit to noise image makes (Stirling) approximate formula, and poisson noise is converted into white Gaussian noise;

Si Di makes (Stirling) approximate formula think, when k is big, Try to achieve by formula (4) is approximate:

$k!\≈\sqrt{2\mathrm{\πk}}{k}^k{e}^{-k}---\left(4\right)$

At this moment, noise is converted into Gaussian distribution from Poisson distribution, and then X ray noise profile expression formula is:

$p\left(k\right)=\frac{1}{\sqrt{2\mathrm{\πh}}}{e}^{[-\frac{{(h-k)}^2}{2h}]}---\left(5\right)$

Then original image meets the white Gaussian noise model, and its mathematical expectation E is asked by formula (6):

$E{\left|\right|c(x_i,y_i)-c(x_j,y_j)\left|\right|}_{2a}^2=E{\left|\right|r(x_i,y_i)-r(x_j,y_j)\left|\right|}_{2a}^2+{2\mathrm{\σ}}^2---\left(6\right)$

Wherein,

is pixel i; J place neighborhood Gauss's weighted euclidean distance; A＞0 is the Gaussian convolution nuclear of standard, and noise variance is σ, the size of the fuzzy de-noising window of size decision of σ;

2. noise level is estimated

Utilize the piecemeal method that image is divided into many sub-pieces in the inventive method earlier, the method for utilizing parallel filtering utilizes particle swarm optimization algorithm to find out maximum noise level in the sub-piece, noise level estimated result as a whole to every sub-block filtering again;

3. fuzzy noise filtering template window

Two window sizes are set; One is neighborhood of pixels window size A * A, and another is the window size B * B of neighborhood of pixels window hunting zone, and promptly selecting the neighborhood size of pixel in the inside, zone of B * B size is A * A; Carry out non local average (the Adaptive Non Localmeans of self-adaptation; ANL-means made in brief note) algorithm, the window of B * B slides in the zone of A * A size, confirms the contribution weights of regional center pixel grey scale according to the similarity in zone;

2). the non local average filter algorithm of bivariate fuzzy self-adaption

The inventive method is utilized the optimal parameter of particle swarm optimization algorithm from wave filter, seeks to make the minimum weighted value of X-ray scanning characteristics of image vector dimensionality reduction reconstruction error, thereby obtains effective denoising method of X-ray scanning image;

The present invention carries out ANL means Filtering Processing in image subblock; Be configured to the individuality of particle swarm optimization algorithm according to the relative program of sample point and its point of proximity; The individuality of population finds global optimum's speed in the process of optimizing then, confirms the global optimum position, finally obtains the most relevant neighbour's partial reconstruction weight matrix; Thereby remove incoherent redundant information effectively, realize effective de-noising of X ray;

1. the non local average filter of bivariate

If c (x y) is a width of cloth observed image, i.e. X ray detected image, and n is that average is 0, variance is σ ²Additive noise, then (x y) is mapped in the observation space (x, y) two-dimentional bounded domain of definition, (x, y) ∈ R through non local average (Non Local means, NL means made in brief note) algorithm input picture r ², x _i, y _iWith neighborhood point x _j, y _jBetween correlation NL (c (x _i, y _i)) obtain by (6) formula:

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\∫e^{-\frac{(G_a*\left({|c(x_i,y_i)-c(x_j,y_j)|}^2\right)\left(0\right))}{h^2}}c(x_j,y_j)\mathrm{dxdy}---\left(6\right)$

Wherein variable be (x, y),

For standard variance is the Gaussian convolution nuclear of a, h is a filtering parameter, and mainly the noise criteria difference by image determines D (x _i, y _i) be NL means conversion coefficient, it according to coordinate (x, the gray-scale value of relative position y) is obtained:

$D(x_i,y_i)={\∫e}^{\frac{(G_a*\left({\left|c(x_i,y_i)c(x_t,y_t)\right|}^2\right)\left(0\right))}{h^2}}{x}^{\′}\left(i\right)y^{\′}\left(i\right)\mathrm{di}---\left(7\right)$

Be expressed as c (i) for digital X ray scanning image in discrete domain; Its pixel i ∈ (x, y), C={c (i); I ∈ I}; I is a collection of pixels in the image, then the pixel NL means algorithm form that is converted into formula (8) according to the most relevant neighbour's partial reconstruction weight matrix ask between pixel i and the neighborhood correlation NL (c (i)) with relevant neighbour's partial reconstruction weight w (i, j) and the gray-scale value c (j) of its neighbor point represent:

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\∈I}{\mathrm{\Σ}}w(i,j)c\left(j\right)---\left(8\right)$

Wherein, in the formula (2), point (x _i, y _i) gray-scale value located is c (i), (x y) changes correction factor γ and transforms its neighbor point (x by c _j, y _j) gray-scale value be c (j), I is that the pixel of X ray detected image is always counted, X is that (x, width value y), Y are c (x, height value y) to c;

Wherein, i={1,2...i...I}={ (1,1), (1,2) ... (x _i, y _i) .... (X, Y) }, at this moment, relevant neighbour's partial reconstruction weight w (i j) is expressed as:

$w(i,j)=\frac{1}{Z\left(i\right)}{e}^{{\left|\right|c(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(9\right)$

$Z\left(i\right)=e^{-{\left|\right|c(\mathrm{NL}(x_{i,},y_i)-c\left(\mathrm{NL}(x_s,y_s)\right)\left|\right|}_{2a}^2/h^2}---\left(10\right)$

Wherein, relevant neighbour's partial reconstruction weight w between pixel i and the neighbor point j (i s) sets up the relation of formula (9), and Z (i) is a weight weighting transformation ratio,

Be the Gauss's weighted euclidean distance based on gray level of two neighborhoods at an i and j place, (i s) is mainly determined by parameter a and h weighted value w; A is the standard deviation of the Gaussian convolution nuclear of standard, and h is the filtering parameter of bivariate NL means, and a is big more; Weights are just more little, show pixel x _i, y _iThe distance center pixel is far away more, and the size of a is by the window size A * A decision of selected neighborhood of pixels; Along with the increase of filtering parameter h, artificial pseudo-shadow weakens gradually, but spatial resolution also can descend, and filtering parameter h is by the standard deviation decision of noise; 2. particle group optimizing

Being located at the local optimum parameter that obtains in the p sub-module is a _pAnd h _p, be the local optimum solution vector, utilize the population optimal algorithm, each particle can be estimated the adaptive value of self-position through certain rule; Each particle can be remembered the own current local desired positions z that finds _pPromptly $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right],$ The local location of particle is z _Pq, z _p∈ z _Pq={ z _P1, z _P2...., z _Pq..z _PQ, z _PqBe particle position in the local group kind, Q is the population that comprises in the population, remembers the desired positions that all individualities find in the colony in addition, is called the z of global optimum _g, promptly $\left[\begin{array}{c}a\\ h\end{array}\right],$ Select in the population best $\left[\begin{array}{c}a\\ h\end{array}\right]$ Solution vector;

Particle swarm optimization algorithm points out that all particles all have speed decision their direction and a distance, is called local optimum speed, uses v _pExpression; Particles are followed current optimal particle and in solution space, are searched for, and in each iteration, particle upgrades position and speed according to following formula:

z _p+1＝z _p+v _p+1 (11)

v _p+1＝w ₁v _p+g ₁(s _pq-z _p)+g ₂(s _gq-z _p) (12)

In the formula, Q is the dimension in target search space, and P representes quantity individual in the population, calculates s _pThe current location adaptive value, v _pBe the flying speed of particle p, s _p∈ (s _P1, s _P2, s _Pq... s _PQ) optimal location that searches up to now for particle, s _PqBe the position of particle process in the population, s _pFor speed is v _pThe time optimal location that searches; s _g=(s _G1+ s _G1+ ... s _Gq+ ... s _GQ) optimal location that searches up to now for whole population and, s _GqBe the optimum position of population current record, g ₁And g ₂Be the self study factor, w ₁Be inertia weight;

3. relevance grade evaluation function

The relevance grade evaluation function is to weigh individual good and bad sign, and its effect is in piecemeal subgraph module, to find out whole optimal parameter a and h; According to the uncertainty of the industrial X-ray scan image noise model singularity as individual model, the relevance grade evaluation function that is adopted is:

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\′}\left|\right|}^2---\left(13\right)$

In the formula, u _PqThe q individuals in p generation in the expression population, D representes quantity individual in the population, the concrete steps of this algorithm optimization are following:

1. initialization population: select threshold epsilon and maximum iteration time u _Max, primary iteration number of times u=0, the particle position scope is z _Min～z _Max, the particle rapidity scope is v _Max～v _Min

2. the individuality in the population is measured: the adaptive value u that measures each particle _Pq, choose this colony's optimal-adaptive value and historical colony optimal-adaptive value compares, obtain the optimal-adaptive value u of colony till current _p, and u _p∈ { u _Pq| u _P1, u _P2... u _Pq... u _PQ, and the value z of access relevant position _g∈ { z _Pq| z _P1, z _P2... z _Pq... z _PQ; Each particle obtains the optimal-adaptive value u till own current through self adaptive value size more in history _g∈ { u _Gq| u _G1, u _G2... u _Gq... u _GQ, and access relevant position z _g∈ { z _Gq| z _G1, z _G2... z _Gq... z _GQ;

3. estimate population U _p, and keep optimum solution;

4. population operation: if U _p＞ε, order is carried out, otherwise end loop;

5. according to rule, upgrade v _p, z _p

6. change evolutionary generation: evolutionary generation adds 1, like the maximum evolutionary generation P of unmet still _MaxAlgorithm goes to this joint step and 2. proceeds.

Images