CN102609904B - 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 development of industrial X-ray inspection technique, the quality of X-ray scanning image has also been proposed to increasing requirement, this just requires effectively to eliminate the noise information producing in real-time testing process.Because X-ray scanning gradation of image interval is narrow, Defect Edge is fuzzy, picture noise is many, defect characteristic covered feature sometimes, affected the effect of detected workpiece being analyzed and being evaluated according to radioscopic image.Along with the continuous renewal 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 more pixel in image, and the amount of redundancy between its pixel increases greatly.Therefore, by fuzzy self-adaption, optimize incoherent amount in operator removal information, and to keep the unchangeability of its inner structure be non-local average de-noising gordian technique urgently to be resolved hurrily.

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 to compare with the gray-scale value of single pixel in image, but the distribution situation of this pixel whole gray scale is around compared, according to the similarity of whole intensity profile, estimate weights, its essence is that the relevance map between original image pixel is processed to image space.It is de-noising service that Buades points out to make full use of redundant information, theoretical analysis and experimental result at Buades research in recent years show, the non local average filter denoising algorithm of tradition is all better than common Image denoising algorithm in subjective and objective performance, as gaussian filtering, anisotropic filtering, total error minimizes, Neighborhood Filtering etc., it originates from Neighborhood Filtering algorithm, that the one of Neighborhood Filtering algorithm is promoted, its weights obtain according to the pixel similarity that around whole area grayscale distributes, when reducing picture noise, there is the ability of very strong maintenance image spatial resolution.While utilizing traditional non local average filter algorithm process complicated image, its calculated amount is larger, and processing speed is slow, and especially, when processing larger image, this problem is more outstanding; In addition, the method can be introduced artificial artifact 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, take the concept of particle group optimizing operator and theory as basis, when having multiple optimization filtering parameters to select, utilize fuzzy control rule to carry out the selection of optimum filtering parameter, by selecting Optimal Window, opened evolutionary search, obtain better treatment effect, faster the ability of speed of convergence and global optimizing.But up to the present, also nobody by this algorithm application in non local average optimization aspect.

For the problems referred to above, the present invention improves traditional NL-means method, and the intrinsic property of algorithm between data found the weighted value that makes reconstruction error minimum, reaches the object that makes error function minimum.Meanwhile, the essence of intelligent optimization algorithm has determined definite distance not relying between data point of weight, effectively reduces the correlativity between Image neighborhood pixel, computational complexity reduces, saved computing time, improved algorithm travelling speed, to adapt to the detection needs in commercial production task.

Summary of the invention

The present invention is studied the de-noising situation of the random noise model producing in industrial X-ray scan image, and the bivariate non local average filter quick noise-eliminating method of employing based on fuzzy self-adaption parameter adjustment realized effective denoising Processing of industrial X-ray scan image.

Fuzzy self-adaption parameter adjustment is effectively to the parallel processing of X ray block image, reduced the correlativity between image block, realize the quick denoising Processing of X-ray scanning image, simultaneously, in invention, introduce x, y bivariate has guaranteed the location invariance of Image Denoising by Use process, and, because adopting particle group optimizing method, it also can obtain good convergence.Thereby the present invention utilizes the intrinsic property of particle swarm optimization algorithm between data, find and in X-ray scanning image, make the weighted value of reconstruction error minimum, thereby make industrial X-ray de-noising image obtain processing speed faster under the condition that guarantees precision.

Below the inventive method is described further, particular content is as follows:

The non local average filter radioscopic image of bivariate noise-eliminating method, feature 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 hypotheses: sampled data place local space is linear, and each sampled point and its Neighbor Points have similarity relation, by 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, according to the similarity of intensity profile, the weight in each neighborhood is set, the map pane that is bumped into of hypothesis, under part is linear condition, minimizes uncorrelated pixel, reconstruct original image;

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

c(x,y)=r(x,y)+n(x,y) (1)

Now need to look for a template to determine the correlativity between c (x, y) and each element of r (x, y), closely and effectively eliminate n (x, y);

Pixel x in c (x, y) image, the c for gray-scale value (i) that y is corresponding represents have

$\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, i is picture point x, gray-scale pixels in y, i ∈ (1,2....I), x ∈ (1,2....X), y ∈ (1,2.....Y) c (i) is the gray-scale value that in X-ray scanning image, i is ordered, c (x, y) reflection image mid point x, the half-tone information at y place, γ is gradation of image modifying factor, therefore, by regulating every bit γ value, regulate the local contrast of X-ray scanning image, make the gradation of image distribution after denoising suitable;

1. random noise model conversion

The main cause of x-ray imaging system image deterioration is system random noise, on X-ray scanning image, show as level or vertical striped, research shows, the generation of X ray and the noise producing with the interaction of material, its distribution thinks that state is discrete, the Markovian process of Time Continuous, all meets Poisson stochastic process in time with on space, if process with independent increments Δ c (t), the probability distribution of its increment is obeyed Poisson distribution, 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 larger, according to very inconvenient in the very complicated practical application of Poisson distribution formula calculation error, therefore, adopt herein and noise image is carried out to 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 larger, By formula (4) is approximate, try to achieve:

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

Now, noise is converted into Gaussian distribution from Poisson distribution, and 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)$

Original image meets white Gaussian noise model, and its mathematical expectation E is required by formula (6):

$E{\left|\right|c\left(x_{i,}{y}_i\right)-c\left(x_{j,}{y}_j\right)\left|\right|}_{2a}^2=E{\left|\right|r\left(x_{i,}{y}_i\right)-r\left(x_{j,}{y}_j\right)\left|\right|}_{2a}^2+2{\mathrm{\σ}}^2---\left(6\right)$

Wherein,

for pixel i, j place neighborhood Gauss's weighted euclidean distance, the Gaussian convolution core that a>0 is standard, noise variance is σ, the size of σ determines the size of fuzzy de-noising window;

2. noise level is estimated

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

Gaussian noise method of estimation has a variety of, and conventional has based on image block and two kinds of methods based on wave filter.Method based on image block is divided into many fritters by image, calculates its each auto-variance, using the average of some minimum variances as estimated result; Method based on wave filter is first carried out once level and smooth to image, then calculates the difference of former figure and level and smooth rear image, error image estimating noise level thus.In the present invention, first utilize Method of Partitioning that image is divided into many sub-blocks, utilize the method for parallel filtering to each sub-block filtering, recycling particle swarm optimization algorithm is found out maximum noise level in sub-block, as a whole noise level estimated result;

The theoretical analysis of Buades research in recent years and experimental result point out that NL-Means algorithm is all better than common Image denoising algorithm in subjective and objective performance, as gaussian filtering, anisotropic filtering, total error minimize, Neighborhood Filtering etc., its weights obtain according to the pixel similarity that around whole area grayscale distributes, and have the ability of very strong maintenance image spatial resolution when reducing picture noise;

3. fuzzy noise Filtering Template window

In order to improve counting yield, select two window concurrent operations, by 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, in the inside, region of B × B size, selecting the Size of Neighborhood of pixel is A × A, carry out non local average (the Adaptive Non Local means of self-adaptation, ANL-means made in brief note) algorithm, the window of B × B slides in the region of A × A size, according to the contribution weights of the similarity definite area center pixel gray scale in region.

2). the non local average filter algorithm of Two-variable Fuzzy self-adaptation

In FANL means algorithm, most crucial problem is to utilize gaussian additive noise level to obtain to make the relevant partial reconstruction weight matrix of reconstruction error minimum, but, this algorithm is to operate for topography's piece, researcher adopts the variable relevant with Euclidean distance to define this weight matrix, by the distance of distance, weigh the size of impact, this makes this algorithm very sensitive to the noise in sample, and this algorithm the convergence speed is fast not in addition.The non local average filter algorithm of Two-variable Fuzzy self-adaptation is based upon fuzzy de-noising window and selects on basis, each little image block is carried out to non local average filter denoising, make to obtain in each module a pair of optimization filtering parameter, and utilize particle swarm optimization algorithm to realize the renewal operation of optimized parameter, thereby realized the Optimization Solution of target.So the present invention utilizes the optimal parameter of particle swarm optimization algorithm from wave filter, find the weighted value that makes X-ray scanning image feature vector dimensionality reduction reconstruction error minimum, thereby obtain effective denoising method of X-ray scanning image.

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

1. the non local average filter of bivariate

If c (x, y) is a width observed image, i.e. X ray detected image, n is that average is 0, variance is σ ²additive noise, input picture r (x, y) is mapped in observation space by non local average (Non Local means, NL means made in brief note) algorithm, (x, y) defines a two-dimentional bounded domain, (x, y) ∈ R ², x _i, y _iwith neighborhood point x _j, y _jbetween correlation NL (c (x _i, y _i)) by (7) formula, obtained:

$\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(7\right)$

Wherein variable is (x, y),

for the Gaussian convolution core that standard variance is a, h is filtering parameter, the main poor decision of noise criteria by image, D (x _i, y _i) be NL means conversion coefficient, it is obtained according to the gray-scale value of the relative position of coordinate (x, y):

$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(8\right)$

For Digitized X-ray scan image, in discrete domain, be expressed as c (i), its pixel i ∈ (x, y), C={c (i), i ∈ I}, I is pixel set in image, correlation NL (c (i)) relevant neighbour's partial reconstruction weight w (i, j) that the form that pixel NL means algorithm is converted into formula (9) according to the most relevant neighbour's partial reconstruction weight matrix is asked between pixel i and neighborhood represents with the gray-scale value c (j) of its neighbor point:

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

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

Wherein, i={1,2...i...I}={ (1,1), (1,2) ... (x _i, y _i) .... (X, Y) }, now, 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\left(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\right||}_{2a}^2/h^2}---\left(10\right)$

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

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

for the Gauss's weighted euclidean distance based on gray level of two neighborhoods at an i and j place, weighted value w (i, s) is mainly determined by parameter a and h, a is the standard deviation of the Gaussian convolution core of standard, and h is the filtering parameter of bivariate NL means, and a is larger, weights are just less, show pixel x _i, y _idistance center pixel is far away, and the size of a is determined by the window size A × A that selectes neighborhood of pixels; Along with the increase of filtering parameter h, artificial artifact weakens gradually, but spatial resolution also can decline, and filtering parameter h is determined by the standard deviation of noise;

2. particle group optimizing

Being located at the local optimum parameter obtaining in p submodule is a _pand h _p, be local optimum solution vector, utilize population optimal algorithm, each particle can be estimated by certain rule the adaptive value of self-position; Each particle can be remembered own current found local desired positions z _p $\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 _pqfor the position of particle in local group kind, Q is the population comprising in population, remembers in addition the desired positions that in colony, all individualities find, and is called the z of global optimum _g, $\left[\begin{array}{c}a\\ h\end{array}\right],$ Select in population best $\left[\begin{array}{c}a\\ h\end{array}\right]$ Solution vector;

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

z _p+1=z _p+v _p+1 (12)

v _p+1=w ₁v _p+g ₁(s _pq-z _p)+g ₂(s _gq-z _p) (13)

In formula, Q is the dimension in target search space, and P represents quantity individual in population, calculates s _pcurrent location adaptive value, v _pfor 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 _pqfor the position of particle process in population, s _pfor speed is v _ptime the 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 _gqfor the optimum position of population current record, g ₁and g ₂for the self study factor, w ₁for inertia weight;

3. relevance grade evaluation function

Relevance grade evaluation function is to weigh individual good and bad sign, and its effect is to find out overall optimal parameter a and h in piecemeal subgraph module.The present invention is the singularity as individual body Model according to the uncertainty of industrial X-ray scan image noise model, and the relevance grade evaluation function adopting 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(14\right)$

In formula, u _pqq the individuality that represents p generation in population, D represents quantity individual in population, the concrete steps of this algorithm optimization are as follows:

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

2. the individuality in 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 by 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. evaluate population U _p, and keeping optimization;

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, as met not yet maximum evolutionary generation P _maxalgorithm goes to this section step and 2. proceeds.

Beneficial effect of the present invention is, in order better to eliminate the impact of the unknown quantum noise existing in industrial X-ray scan image, proposed to transfer not tractable quantum noise model to common gaussian additive noise model, use the size of fuzzy operation selective filter window, find the radioscopic image noise-eliminating method of the Two-variable Fuzzy self-adaptation nonlinear average filter of the relevant weight matrix that makes error function minimum.In the present invention, introduce particle group optimizing filtering parameter, and then partial reconstruction weight matrix, reduced the impact of local correlations on sample data, improved algorithm the convergence speed, the speed and the precision that have improved the denoising of industrial X-ray scan image, be applicable to the processing to the uncertain X-ray scanning image of noise model.

Accompanying drawing explanation

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

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

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

Embodiment

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

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

Non local average filter algorithm has a hypotheses: sampled data place local space is linear, and each sampled point and its Neighbor Points have similarity relation, by 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, according to the similarity of intensity profile, the weight in each neighborhood is set, the map pane that is bumped into of hypothesis, under part is linear condition, minimizes uncorrelated pixel, reconstruct original image;

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

c(x,y)=r(x,y)+n(x,y) (1)

Now need to look for a template to determine the correlativity between c (x, y) and each element of r (x, y), closely and effectively eliminate n (x, y);

Pixel x in c (x, y) image, the c for gray-scale value (i) that y is corresponding represents have

$\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, i is picture point x, gray-scale pixels in y, i ∈ (1,2....I), x ∈ (1,2....X), y ∈ (1,2.....Y) c (i) is the gray-scale value that in X-ray scanning image, i is ordered, c (x, y) reflection image mid point x, the half-tone information at y place, γ is gradation of image modifying factor, therefore, by regulating every bit γ value, regulate the local contrast of X-ray scanning image, make the gradation of image distribution after denoising suitable;

1. random noise model conversion

The main cause of x-ray imaging system image deterioration is system random noise, on X-ray scanning image, show as level or vertical striped, research shows, the generation of X ray and the noise producing with the interaction of material, its distribution thinks that state is discrete, the Markovian process of Time Continuous, all meets Poisson stochastic process in time with on space, if process with independent increments Δ c (t), the probability distribution of its increment is obeyed Poisson distribution, 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 larger, according to very inconvenient in the very complicated practical application of Poisson distribution formula calculation error, therefore, adopt herein and noise image is carried out to 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 larger, By formula (4) is approximate, try to achieve:

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

Now, noise is converted into Gaussian distribution from Poisson distribution, and 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)$

Original image meets white Gaussian noise model, and its mathematical expectation E is required by formula (6):

$E{\left|\right|c\left(x_{i,}{y}_i\right)-c\left(x_{j,}{y}_j\right)\left|\right|}_{2a}^2=E{\left|\right|r\left(x_{i,}{y}_i\right)-r\left(x_{j,}{y}_j\right)\left|\right|}_{2a}^2+2{\mathrm{\σ}}^2---\left(6\right)$

Wherein,

for pixel i, j place neighborhood Gauss's weighted euclidean distance, the Gaussian convolution core that a>0 is standard, noise variance is σ, the size of σ determines the size of fuzzy de-noising window;

2. noise level is estimated

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

Gaussian noise method of estimation has a variety of, and conventional has based on image block and two kinds of methods based on wave filter.Method based on image block is divided into many fritters by image, calculates its each auto-variance, using the average of some minimum variances as estimated result; Method based on wave filter is first carried out once level and smooth to image, then calculates the difference of former figure and level and smooth rear image, error image estimating noise level thus.In the present invention, first utilize Method of Partitioning that image is divided into many sub-blocks, utilize the method for parallel filtering to each sub-block filtering, recycling particle swarm optimization algorithm is found out maximum noise level in sub-block, as a whole noise level estimated result;

The theoretical analysis of Buades research in recent years and experimental result point out that NL-Means algorithm is all better than common Image denoising algorithm in subjective and objective performance, as gaussian filtering, anisotropic filtering, total error minimize, Neighborhood Filtering etc., its weights obtain according to the pixel similarity that around whole area grayscale distributes, and have the ability of very strong maintenance image spatial resolution when reducing picture noise;

3. fuzzy noise Filtering Template window

In order to improve counting yield, select two window concurrent operations, by 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, in the inside, region of B × B size, selecting the Size of Neighborhood of pixel is A × A, carry out non local average (the Adaptive Non Local means of self-adaptation, ANL-means made in brief note) algorithm, the window of B × B slides in the region of A × A size, according to the contribution weights of the similarity definite area center pixel gray scale in region.

2). the non local average filter algorithm of Two-variable Fuzzy self-adaptation

In FANL means algorithm, most crucial problem is to utilize gaussian additive noise level to obtain to make the relevant partial reconstruction weight matrix of reconstruction error minimum, but, this algorithm is to operate for topography's piece, researcher adopts the variable relevant with Euclidean distance to define this weight matrix, by the distance of distance, weigh the size of impact, this makes this algorithm very sensitive to the noise in sample, and this algorithm the convergence speed is fast not in addition.The non local average filter algorithm of Two-variable Fuzzy self-adaptation is based upon fuzzy de-noising window and selects on basis, each little image block is carried out to non local average filter denoising, make to obtain in each module a pair of optimization filtering parameter, and utilize particle swarm optimization algorithm to realize the renewal operation of optimized parameter, thereby realized the Optimization Solution of target.So the present invention utilizes the optimal parameter of particle swarm optimization algorithm from wave filter, find the weighted value that makes X-ray scanning image feature vector dimensionality reduction reconstruction error minimum, thereby obtain effective denoising method of X-ray scanning image.

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

1. the non local average filter of bivariate

If c (x, y) is a width observed image, i.e. X ray detected image, n is that average is 0, variance is σ ²additive noise, input picture r (x, y) is mapped in observation space by non local average (Non Local means, NL means made in brief note) algorithm, (x, y) defines a two-dimentional bounded domain, (x, y) ∈ R ², x _i, y _iwith neighborhood point x _j, y _jbetween correlation NL (c (x _i, y _i)) by (7) formula, obtained:

$\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(7\right)$

Wherein variable is (x, y),

for the Gaussian convolution core that standard variance is a, h is filtering parameter, the main poor decision of noise criteria by image, D (x _i, y _i) be NL means conversion coefficient, it is obtained according to the gray-scale value of the relative position of coordinate (x, y):

$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(8\right)$

For Digitized X-ray scan image, in discrete domain, be expressed as c (i), its pixel i ∈ (x, y), C={c (i), i ∈ I}, I is pixel set in image, correlation NL (c (i)) relevant neighbour's partial reconstruction weight w (i, j) that the form that pixel NL means algorithm is converted into formula (8) according to the most relevant neighbour's partial reconstruction weight matrix is asked between pixel i and neighborhood represents with the gray-scale value c (j) of its neighbor point:

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

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

Wherein, i={1,2...i...I}={ (1,1), (1,2) ... (x _i, y _i) .... (X, Y) }, now, 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\left(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\right||}_{2a}^2/h^2}---\left(10\right)$

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

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

for the Gauss's weighted euclidean distance based on gray level of two neighborhoods at an i and j place, weighted value w (i, s) is mainly determined by parameter a and h, a is the standard deviation of the Gaussian convolution core of standard, and h is the filtering parameter of bivariate NL means, and a is larger, weights are just less, show pixel x _i, y _idistance center pixel is far away, and the size of a is determined by the window size A × A that selectes neighborhood of pixels; Along with the increase of filtering parameter h, artificial artifact weakens gradually, but spatial resolution also can decline, and filtering parameter h is determined by the standard deviation of noise;

2. particle group optimizing

Being located at the local optimum parameter obtaining in p submodule is a _pand h _p, be local optimum solution vector, utilize population optimal algorithm, each particle can be estimated by certain rule the adaptive value of self-position; Each particle can be remembered own current found local desired positions z _p $\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 _pqfor the position of particle in local group kind, Q is the population comprising in population, remembers in addition the desired positions that in colony, all individualities find, and is called the z of global optimum _g, $\left[\begin{array}{c}a\\ h\end{array}\right],$ Select in population best $\left[\begin{array}{c}a\\ h\end{array}\right]$ Solution vector;

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

z _p+1=z _p+v _p+1 (12)

v _p+1=w ₁v _p+g ₁(s _pq-z _p)+g ₂(s _gq-z _p) (13)

In formula, Q is the dimension in target search space, and P represents quantity individual in population, calculates s _pcurrent location adaptive value, v _pfor 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 _pqfor the position of particle process in population, s _pfor speed is v _ptime the 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 _gqfor the optimum position of population current record, g ₁and g ₂for the self study factor, w ₁for inertia weight;

3. relevance grade evaluation function

Relevance grade evaluation function is to weigh individual good and bad sign, and its effect is to find out overall optimal parameter a and h in piecemeal subgraph module.The present invention is the singularity as individual body Model according to the uncertainty of industrial X-ray scan image noise model, and the relevance grade evaluation function adopting 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(14\right)$

In formula, u _pqq the individuality that represents p generation in population, D represents quantity individual in population, the concrete steps of this algorithm optimization are as follows:

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

2. the individuality in 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 by 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. evaluate population U _p, and keeping optimization;

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, as met not yet maximum evolutionary generation P _maxalgorithm goes to this section step and 2. proceeds.

The present invention is directed to the features such as in electric field commercial unit Non-Destructive Testing process, X-ray scanning gradation of image interval is narrow, Defect Edge is fuzzy, picture noise is many, defect characteristic is submerged sometimes and carry out Denoising Study, respectively with horizontal quantum noise, vertical quantum noise, random noise, the noise of crossing enhancing processing generation is that example is studied, the scan image of the multiple format image arbitrary sizes such as input BMP, JPG, PNG.First by user, according to denoising index request, auto-selecting parameter (target image) is proposed, as user does not provide by system automatic discrimination, and the quantum noise of the image of input is converted to Gaussian noise, then, high phase noise is carried out to grade classification, according to the distribution of noise model level, utilize fuzzy control rule that different windows is set, and different filtering parameters is set in each subwindow, utilize Particle Swarm Optimization to carry out the selection of global optimum's filtering parameter, closely view picture X-ray scanning image is carried out to FANL denoise algorithm, and with least mean-square error (Minimum Mean Squared Error, be MMSE), Y-PSNR (peak signal noise ratio, be PSNR) and denoising time t used evaluation system denoising performance.As shown in Figure 2, specific implementation method is as follows for system architecture diagram:

1. noise level modeling and template window size Selection

The main cause of x-ray imaging system image deterioration is system random noise.Research shows, the generation of X ray and with the interaction of material, all meet Poisson stochastic process in time with on space.For rapid X-ray imaging system, because the time shutter is short, the quantum noise that X ray produces is more outstanding, has had a strong impact on the quality of image.Therefore,, to the X-ray scanning image of input, due to its noise model the unknown, therefore, on original image basis, add average μ ₀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, now, the image function model that utilizes this base of a fruit to make (Stirling) formula (4) add after the even noise of Gauss X-ray scanning is similar to, poisson noise is converted into white Gaussian noise, the variances sigma of obtaining again the overall white Gaussian noise of the rear image of conversion, σ value affects system performance.

According to fuzzy noise Filtering Template window selection principle, need to add image after Gauss's noise processed to the radioscopic image after detecting carries out fuzzy block parallel and calculates each local optimum parameter, therefore, two windows need to be set, neighborhood of pixels window size A × A, window B × the B of neighborhood window hunting zone, the image that common noise is larger, generally gets A=7, B=23, for low noise image, get A=3, B=7 just can meet noise reduction substantially.Therefore, noise level arranges the template window size of non local average de-noising, and for the noise of input, first adding parameter is (u ₀, σ ₀) white Gaussian noise after conversion noise model, according to the difference of the high phase noise level σ in image function after conversion, be divided into five grades, be low noise σ <0.2, compared with low noise 0.2< σ <0.5, medium noise 0.5< σ <0.8, higher noise 0.8< σ <1 and five domain intervals of strong noise σ >1, for different noises, different processing template windows is set, in each template window, adjust adaptively filtering parameter, under varying level, the class level of obscures window as shown in Figure 3.

2. in template window, local optimum parameter is selected

(1) local parameter a

In non local average filter, the weighted value w (i, j) of pixel i and vicinity points j is mainly determined by parameter a and h, the Gaussian convolution core that a is standard, h is the filtering parameter of bivariate NL means, a is larger for Gaussian convolution core, and weights are just less, show pixel x _i, y _idistance center pixel is far away, and the size of a is determined by the window size A × A that selectes neighborhood of pixels; The level and smooth degree of the value of a and image itself also has certain relation, when a value is larger, comprised more noise spot, but also can will in image, should comprise to come in by the point different from current pixel value simultaneously, therefore a value is not the bigger the better, and a often gets in the scope of 1～10 times of neighborhood window size; Take into account the impact of different noise levels simultaneously, get 10～15 times of the poor σ of noise criteria as with reference to value, consider that all and current pixel weighting absolute difference is no more than the pixel of a σ.Therefore,, when the optimized parameter of selecting is discontented while being enough to upper condition, give up the local optimum parameter value of having selected.

(2) local parameter h

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

When adding while making an uproar after processing that the distribution of pixel meets Gaussian distribution in image, its noise weight function w (i, j) also meeting average is 0, standard deviation is the typical Gaussian function of h/2, therefore, by changing h value, regulate weight w (i, j) distribution, w (i, j) meet the feature of Gaussian distribution, 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 (σ) having superposeed, these differences have determined i, the Gauss weighted euclidean distance value of j point-to-point transmission based on gray level, within having a certain proportion of pixel to drop on to be no more than the poor β of noise criteria scope doubly with the weighting absolute difference of current considered pixel value.Setting is ignored weights and is ignored noise pixel number percent with calculating during Weighted distance, and all weights sums that are greater than a certain multiple of standard deviation when weights distance estimate, set the weights ignored and, ask and obtain this multiple value.It is β that the weighting absolute difference of now supposing pixel value is no more than the poor multiple of noise criteria, and not uncared-for distance should meet:

${\left|\right|c\left({\mathrm{NL}}_i\right)-c\left({\mathrm{NL}}_j\right)\left|\right|}_{2a}^2\&le;\frac{\mathrm{\&beta;h}}{\sqrt{2}}---\left(15\right)$

The approximation relation that obtains filtering parameter and noise level is:

$h=\sqrt{2}{\mathrm{\&beta;\&sigma;}}^2---\left(16\right)$

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

$h=[\sqrt{2}\mathrm{\&beta;}+\frac{\max \left(O({\mathrm{\&sigma;}}_c-{\mathrm{\&sigma;}}_0)\right)}{10}]{\mathrm{\&sigma;}}^2---\left(17\right)$

σ _cfor the average of detected image poor, σ ₀for normative reference difference limen value, when the parameter in employing formula (17) 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; H parameter value is excessive, can smooth out the image difference that some exceed noise scope, thereby causes the detailed information such as image is too level and smooth, loss edge.For gray scale, be the X-ray scanning image of 256 grades, be less than σ ₀standard deviation less on h value impact, and when noise level is higher than σ ₀in time, affect larger on h.

3. the filtering parameter solution vector based on particle group optimizing is selected

When radioscopic image is carried out to denoising Processing, first will carry out piecemeal to entire image, recycling ANL means algorithm is obtained the local optimum parametric solution vector of each sub-block.Because X-ray production apparatus scans the noise model the unknown in obtaining, therefore different to the number of different x-ray image block, local optimum parameter in each image subblock is also not quite similar, and therefore, adopts particle swarm optimization algorithm in X-ray scanning image, to find global optimum's parametric solution vector.

In the X-ray scanning image subblock of the present invention after Gaussian noise model conversion, according to the correlativity of each pixel gray-scale value, carry out characteristic information extraction, wherein, neighborhood window size A × A, the window size B × B of hunting zone, now, not only between consecutive point, be correlated with, non-conterminously in setting range, also exist certain correlativity, weighted value w (i, j) is according to pixel i, correlativity between j obtains, and shows the local optimum parameter a that obtains each sub-block in each sub-block _i, h _i, recycling particle swarm optimization algorithm is obtained overall optimal parameter a, h.

If input radioscopic image is divided into P sub-block, consider the correlativity between pixel, between each sub-block, also exist certain correlativity, 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 sample point, by distance-taxis, k and c before selecting _pnearest point is as its neighbor point.

By the Neighbor Points of each sample point, calculated the partial reconstruction weight matrix of this sample point $\left[\begin{array}{c}{a}_p\\ h_p\end{array}\right].$ With the Neighbor Points of each proper vector, this proper vector is rebuild, ask for the neighbour's partial reconstruction weight matrix that makes reconstruction error minimum.

It is as follows that error of fitting function is rebuild in definition:

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

Wherein, c _q(q=1,2 ..., k) be c _pq in a sample Neighbor Points, u _pqc _pwith c _qbetween weights.

When meeting following two constraint condition, by minimum error function, obtain partial reconstruction weight matrix, by the Neighbor Points of sample point, construct optimum W matrix and make error function value reach minimum.

1) each data point c _qall can only be represented by its reference point, if c _qnot reference point, u _pq=0;

2) every a line of weight matrix and be 1, meet normalization constraint

Asking in the process of optimum U matrix, U rebuilds the error function that error of fitting Function Solution vector forms, primary using the reconstruction weight vector set of sample point and its reference point as particle cluster algorithm in the present invention, then the individuality in population finds global optimum position and optimum velocity in the process of optimizing, finally obtains neighbour's partial reconstruction weight matrix U.Wherein, u _pqq the individuality that represents p generation in population, P represents quantity individual in population.The concrete steps of this algorithm optimization are as follows, and its algorithm flow as shown in Figure 1.

In formula, u _pqq the individuality that represents p generation in population, D represents quantity individual in 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, particle position scope is z _min～z _max, particle rapidity scope is v _max～v _min;

2) individuality in 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 by 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) evaluate population U _p, and keeping optimization;

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, as met not yet maximum evolutionary generation P _maxalgorithm goes to step 2) proceed.

According to above step, ask for optimum neighbour's partial reconstruction weight matrix U, for the weights U (P) that keeps upper step to ask for is 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) speed is upgraded operation: v _p+1=w ₁v _p+ g ₁(s _pq-z _p)+g ₂(s _gq-z _p)

In formula, q ∈ { q ₁, q ₂... Q}, Q is the dimension in target search space, P is the number of data point in sample, calculates s _pcurrent location adaptive value, v _pfor 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 _pqfor the position of particle process in population, s _pfor speed is v _ptime the 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 _gqfor the optimum position of population current record, g ₁and g ₂for the self study factor, w ₁for inertia weight.Therefore, minimizing vector corresponding to minimum Q eigenwert of getting f (u) in error of fitting function is column vector composition matrix [a _p, h _p] ^t, the column vector of U is the vector representation of Q dimension space.

4. similarity measurement

In the present invention, the quantum noise in the industrial X-ray scan image of input is converted to Gaussian noise, closely the imagery exploitation FANL algorithm with Gaussian noise after conversion is carried out to quick denoising.Therefore, the radioscopic image with Gaussian noise after the industrial X-ray detected image with uncertain quantum noise detecting and the noise model that adds low noise level are changed carries out similarity measurement, and the present invention adopts following similar Measurement Method:

(1) mutual information tolerance

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

Suppose that the industrial X-ray detected image with uncertain quantum noise detecting is image C R, the radioscopic image with Gaussian noise after changing is image C, and C is C by selecting the sub collective drawing after neighborhood window piecemeal ₁, C ₂..C _a..C _a, C is divided into C by search neighborhood window ₁, C ₂..C _b.C _b, CR, C is any real number, CR, the marginal probability density of C is respectively P _cRand P (cr) _c(c), their joint probability distribution is expressed as P _aB(a, b), represents two degrees of correlation between signal by simple crosscorrelation information, has

$I\left(\mathrm{CR},C\right)=\underset{\mathrm{cr},c}{\mathrm{\&Sigma;}}{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(18\right)$

Simple crosscorrelation and entropy have close relationship, and formula (18) is write the form of an accepted way of doing sth (19):

I(CR,C)=H(CR,C)=H(CR)+H(C)-H(CR,C)=H(C)H(C/CR) (19)

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

$H\left(\mathrm{CR}\right)=-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{pi}\log \mathrm{pi}---\left(20\right)$

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

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

$H(\mathrm{CR},C)=-\underset{\mathrm{cr},c}{\mathrm{\&Sigma;}}{P}_{\mathrm{CRC}}(\mathrm{cr},c)\log P_{\mathrm{CRC}}(\mathrm{cr},c)---\left(22\right)$

If the gray-scale value that m (e) is former detected image is at a gray-scale value at s place, e is the pixel in image, Γ (T (e)) carries out the gray-scale value after Γ conversion after the dry acoustic model conversion of original image at an e place, gray-scale value size is set between 0～255, 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, the gray scale that is denoted as two width figure compositions is h to the number of (Γ (e), m (T (e))) _α(m), the image C R that two width are associated and the joint distribution probability of C are Γ

marginal distribution probability

with

so cross correlation function is:

$I(\&PartialD;)=\underset{\mathrm{cr},c}{\mathrm{\&Sigma;}}{P}_{\mathrm{CRC},\&PartialD;}(\mathrm{cr},c){\log }_2\frac{P_{\mathrm{CRC},\&PartialD;}(\mathrm{cr},c)}{P_{\mathrm{CR},\&PartialD;}\left(\mathrm{cr}\right)P_C\left(c\right)}---\left(23\right)$

Optimum mutual information parameter

have

${\&PartialD;}^{*}=\mathrm{arc}\max I\left(\mathrm{\&alpha;}\right)---\left(24\right)$

Obtain respectively neighborhood window C ₁, C ₂..C _a..C _awith search window C ₁, C ₂..C _b.C _boptimum mutual information parameter

with

optimal value ${{\&PartialD;}_o}^{*}=\max \{{\&PartialD;}_a^{*},{\&PartialD;}_b^{*}\}.$

(2) cosine distance metric

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

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

$l=\underset{w=1}{\overset{N}{\mathrm{\&Sigma;}}}\min (F_w,Q_w)---\left(25\right)$

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

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

Trying to achieve cosine distance is:

l(F,Q)=F ^T*Q/(||F||*||Q||) (27)

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

For 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{\&Sigma;}}}*l_w(F,Q)---\left(28\right)$

Expression is formed by m Feature Combination, wherein u _wrepresent the weight coefficient of w feature, it represents the importance of w feature, generally gets each u _wequate S _w(F, Q) represents the similarity measurement functional value of w feature.

5. the performance evaluation of algorithm

The present invention adopts conventional least mean-square error MMSE(Minimum Means Squared Error in denoising performance evaluation), Y-PSNR PSNR and response time t carry out the quality of objective evaluation denoising performance.

The intensity of variation of view data after MMSE evaluation denoising, the value of MMSE is less, illustrates that denoising effect is more accurate.The noise model of former X ray detected image is quantum noise, and the noise square error in the Gaussian noise model drawing after being changed is σ ₀, image is after the method denoising Processing in the present invention, and the information entropy of new images each point pixel after known de-noising, by the Pixel Information entropy comparison of the pixel entropy of image after de-noising and former X-ray scanning image each point, obtains overall meansquaredeviationσ ₁, MMSE=min{ σ ₀, σ ₁, the situation of the average de-noising level of image each point after denoising in MMSE reflection.For a little bigger impact on denoising effect of difference in key diagram picture, in the present invention, introduce PSNR as another objective standard of evaluating image denoising effect quality, for weighing the situation of change between image maximum of points and noise, it is CR that known detection obtains X-ray scanning image, it is C that CR is converted to image through noise model, obtains image R after C de-noising ^*, ideally expect that the not Noise obtaining is R at image, shown in PSNR is defined as follows:

$\mathrm{PSNR}=10\&times;\log \left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(29\right)$

$\mathrm{MMSE}=\min \left\{\right|\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}({\left(R^x\right)}^n-R^n)}{N}\left|\right\}---\left(30\right)$

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

Weighing the response time in denoising system performance index, it is also an index of ignoring, in the present invention owing to adopting each module parallel computation after piecemeal, thereby be better than traditional NL means denoise algorithm in the processing time, but consider precision requirement aspect, in the present invention, utilize x (n), y (n) describes the positional information of pixel n in image, thereby can process more accurately the denoising situation of appointed area.

While applying traditional non local average filter denoising algorithm processing complicated image, its calculated amount is larger, and processing speed is slow, and especially, when processing larger image, this problem is more outstanding; In addition, traditional non local average filter method can be introduced artificial artifact at the smooth region of image, and image thickens, and spatial resolution is affected.The non local average filter denoising algorithm of fuzzy bivariate proposing in the present invention utilizes the mode of piecemeal that complicated image is carried out to piecemeal, again sub-block is carried out to auto-adaptive parameter concurrent operation, when reducing image complexity, saved computing time, reached good denoising performance.

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 hypotheses: sampled data place local space is linear, and each sampled point and its Neighbor Points have similarity relation, by 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 the map pane that is bumped into of hypothesis, under part is linear condition, minimizes uncorrelated pixel, reconstruct original image;

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

c(x,y)=r(x,y)+n(x,y) (1)

Now need to look for a template to determine the correlativity between c (x, y) and each element of r (x, y), closely and effectively eliminate n (x, y);

Pixel x in c (x, y) image, the c for gray-scale value (i) that y is corresponding represents have

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

Wherein, I=X*Y is the size dimension of image, i is picture point x, gray-scale pixels in y, i ∈ (1,2....I), x ∈ (1,2....X), y ∈ (1,2.....Y) c (i) is the gray-scale value that in X-ray scanning image, i is ordered, c (x, y) reflection image mid point x, the half-tone information at y place, γ is gradation of image modifying factor, therefore, by regulating every bit γ value, regulate the local contrast of X-ray scanning image, make the gradation of image distribution after denoising suitable;

1. random noise model conversion

The main cause of x-ray imaging system image deterioration is system random noise, on X-ray scanning image, show as level or vertical striped, the generation of X ray and the noise producing with the interaction of material, its distribution is that state is discrete, the Markovian process of Time Continuous, all meet Poisson stochastic process in time with on space, if process with independent increments Δ c (t), the probability distribution of its increment is obeyed Poisson distribution, has:

$p\left(k\right)\&ap;\frac{{\mathrm{\&lambda;}}^k}{k!}{e}^{-\mathrm{\&lambda;}}---\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 larger, according to very inconvenient in the very complicated practical application of Poisson distribution formula calculation error, 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 larger, By formula (4) is approximate, try to achieve:

$k!\&ap;\sqrt{2\mathrm{\&pi;k}}{k}^k{e}^{-k}---\left(4\right)$

Now, noise is converted into Gaussian distribution from Poisson distribution, and X ray noise profile expression formula is:

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

Original image meets white Gaussian noise model, and its mathematical expectation E is required by formula (6):

$E{\left|\right|c\left(x_{i,}{y}_i\right)-c\left(x_{j,}{y}_j\right)\left|\right|}_{2a}^2=E{\left|\right|r\left(x_{i,}{y}_i\right)-r\left(x_{j,}{y}_j\right)\left|\right|}_{2a}^2+2{\mathrm{\&sigma;}}^2---\left(6\right)$

Wherein,

for pixel i, j place neighborhood Gauss's weighted euclidean distance, the Gaussian convolution core that a>0 is standard, noise variance is σ, the size of σ determines the size of fuzzy de-noising window;

2. noise level is estimated

In the inventive method, first utilize Method of Partitioning that image is divided into many sub-blocks, utilize the method for parallel filtering to each sub-block filtering, recycling particle swarm optimization algorithm is found out maximum noise level in sub-block, as a whole noise level estimated result;

3. fuzzy noise Filtering Template window

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, in the inside, region of B × B size, selecting the Size of Neighborhood of pixel is A × A, carry out the non local average algorithm of self-adaptation: Adaptive Non Local means, ANL-means made in brief note, and the window of B × B slides in the region of A × A size, according to the contribution weights of the similarity definite area center pixel gray scale in region;

2). the non local average filter algorithm of Two-variable Fuzzy self-adaptation

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

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

1. the non local average filter of bivariate

If c (x, y) is a width observed image, i.e. X ray detected image, n is that average is 0, variance is σ ²additive noise, input picture r (x, y) by non local average algorithm: Non Local means, NL means made in brief note, is mapped in observation space, (x, y) defines a two-dimentional bounded domain, (x, y) ∈ R ², x _i, y _iwith neighborhood point x _j, y _jbetween correlation NL (c (x _i, y _i)) by (7) formula, obtained:

$\mathrm{NL}(c(x_i,y_i)=\frac{1}{D(x_i,y_i)}\&Integral;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(7\right)$

Wherein variable is (x, y),

for the Gaussian convolution core that standard variance is a, h is filtering parameter, the main poor decision of noise criteria by image, D (x _i, y _i) be NL means conversion coefficient, it is obtained according to the gray-scale value of the relative position of coordinate (x, y):

$D(x_i,y_i)={\&Integral;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}^{\&prime;}\left(i\right)y^{\&prime;}\left(i\right)\mathrm{di}---\left(8\right)$

For Digitized X-ray scan image, in discrete domain, be expressed as c (i), its pixel i ∈ (x, y), C={c (i), i ∈ I}, I is pixel set in image, correlation NL (c (i)) relevant neighbour's partial reconstruction weight w (i, j) that the form that pixel NL means algorithm is converted into formula (9) according to the most relevant neighbour's partial reconstruction weight matrix is asked between pixel i and neighborhood represents with the gray-scale value c (j) of its neighbor point:

$\mathrm{NL}\left(c\left(i\right)\right)=\underset{s\&Element;I}{\mathrm{\&Sigma;}}w(i,j)c\left(j\right)---\left(9\right)$

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

Wherein, i={1,2...i...I}={ (1,1), (1,2) ... (x _i, y _i) .... (X, Y) }, now, 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\left(\mathrm{NL}(x_{i,},y_i)c\left(\mathrm{NL}(x_s,y_s)\right)\right||}_{2a}^2/h^2}---\left(10\right)$

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

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

for the Gauss's weighted euclidean distance based on gray level of two neighborhoods at an i and j place, weighted value w (i, s) is mainly determined by parameter a and h, a is the standard deviation of the Gaussian convolution core of standard, and h is the filtering parameter of bivariate NL means, and a is larger, weights are just less, show pixel x _i, y _idistance center pixel is far away, and the size of a is determined by the window size A × A that selectes neighborhood of pixels; Along with the increase of filtering parameter h, artificial artifact weakens gradually, but spatial resolution also can decline, and filtering parameter h is determined by the standard deviation of noise;

2. particle group optimizing

Being located at the local optimum parameter obtaining in p submodule is a _pand h _p, be local optimum solution vector, utilize population optimal algorithm, each particle can be estimated by certain rule the adaptive value of self-position; Each particle can be remembered own current found local desired positions z _p $\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 _pqfor the position of particle in local group kind, Q is the population comprising in population, remembers in addition the desired positions that in colony, all individualities find, and is called the z of global optimum _g, $\left[\begin{array}{c}a\\ h\end{array}\right],$ Select in population best $\left[\begin{array}{c}a\\ h\end{array}\right]$ Solution vector;

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

z _p+1=z _p+v _p+1 (12)

v _p+1=w ₁v _p+g ₁(s _pq-z _p)+g ₂(s _gq-z _p) (13)

In formula, Q is the dimension in target search space, and P represents quantity individual in population, calculates s _pcurrent location adaptive value, v _pfor 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 _pqfor the position of particle process in population, s _pfor speed is v _ptime the 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 _gqfor the optimum position of population current record, g ₁and g ₂for the self study factor, w ₁for inertia weight;

3. relevance grade evaluation function

Relevance grade evaluation function is to weigh individual good and bad sign, and its effect is to find out overall optimal parameter a and h in piecemeal subgraph module; Singularity according to the uncertainty of industrial X-ray scan image noise model as individual body Model, the relevance grade evaluation function adopting is:

$f\left(u\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}{\left|\right|v_p-\underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{u}_{\mathrm{pq}}{z_{\mathrm{pq}}}^{\&prime;}\left|\right|}^2---\left(14\right)$

In formula, u _pqq the individuality that represents p generation in population, D represents quantity individual in population, the concrete steps of this algorithm optimization are as follows:

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

2. the individuality in 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 by 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. evaluate population U _p, and keeping optimization;

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, as met not yet maximum evolutionary generation P _maxalgorithm goes to this section step and 2. proceeds.

Images