CN109191399B

CN109191399B - Magnetic resonance image denoising method based on improved multi-path matching pursuit algorithm

Info

Publication number: CN109191399B
Application number: CN201810999368.1A
Authority: CN
Inventors: 范虹; 杨晶; 姚若侠
Original assignee: Shaanxi Normal University
Current assignee: Shaanxi Normal University
Priority date: 2018-08-29
Filing date: 2018-08-29
Publication date: 2022-01-28
Anticipated expiration: 2038-08-29
Also published as: CN109191399A

Abstract

A magnetic resonance image denoising method, involving image processing. This method introduces adaptive genetic algorithm optimization in the process of solving multiple candidate atoms that match local image features in each iteration of the multi-path matching pursuit algorithm, and organically combines the advantages of the adaptive genetic algorithm and the multi-path matching and pursuit algorithm. It can avoid the defect that the genetic algorithm is easy to fall into the local optimum, obtain the best matching parameters with higher precision, and effectively reduce the calculation amount of the multi-path matching pursuit algorithm, and overcome the disadvantage that it cannot be popularized and applied due to the large amount of calculation.

Description

Magnetic resonance image denoising method based on improved multi-path matching pursuit algorithm

Technical Field

The disclosure belongs to the field of image processing, and particularly relates to a magnetic resonance image denoising method.

Background

Magnetic Resonance Imaging (MRI) is a technique for imaging human tissues by using a static Magnetic field and a radio frequency Magnetic field, and the formed image has a vivid anatomical structure, is more clear and intuitive to the lesion part, and has the characteristic of high soft tissue contrast ratio. In recent years, it has been widely used clinically. However, in the scanning process of the MR image, due to the influence of stimulated human tissues and electronic elements, random noise always exists in the image, and a series of problems of uneven gray scale, fuzzy boundary and the like appear, which brings great inconvenience to the processing of subsequent images. Therefore, removing noise in the original MR image has been a challenging problem, and many denoising methods and techniques are applied to denoising the MR image.

The conventional MR image denoising method can be performed in the spatial domain and the frequency domain. Median filtering is a typical nonlinear denoising method based on ordering statistical theory in the spatial domain. The method is characterized in that the pixels in the local area are sorted according to the gray level, the median value in the area is taken as the gray value of the current area, the algorithm principle is simple, and the edge characteristic of the image can be effectively maintained. The main idea of the method is to remove the noise information in the image as a high-frequency part, thereby retaining the low-frequency information in the image to achieve the purpose of denoising. However, there is some detail and edge information in the image that also exists in the high frequency part, so that the processing loses the image information. In recent years, image denoising using the idea of sparse representation has been widely used. In 1993, Mallat and Zhang firstly propose sparse representation idea to acquire sparse expression of an image on a certain transform domain, and propose a Matching Pursuit (MP) algorithm to realize signal sparse decomposition. The denoising principle of the sparse decomposition theory is that atoms are adopted to represent effective information of an image, and noise is random and has no fixed structure, so that the effective information and the noise of a noisy image can be separated, and the purpose of denoising is achieved. In the study of image denoising algorithm based on sparse representation, the core task is the generation of an atom dictionary and the sparse reconstruction algorithm. In 1995, Mallat first applied the MP algorithm to an image to obtain a sparse representation of the image, but since the projections of the residual (signal) in the algorithm on the selected atoms are not necessarily orthogonal, the result of each iteration is not optimal, and multiple iterations are needed to converge. To address this shortcoming, Pati et al orthogonalize each selected atom based on the MP algorithm, and propose an Orthogonal Matching Pursuit (OMP) algorithm, which accelerates the convergence of the algorithm while ensuring the accuracy of the algorithm. Since then, many scholars have proposed some improved algorithms based on the OMP algorithm. The algorithm only has one candidate atomic support set after each iteration, and if an error occurs in a certain iteration process, the final reconstruction result is influenced. Kwon et al were inspired by the problem of combining tree searches, which proposed an algorithm called multipath matching pursuit. The algorithm combines the ideas of a combination method and a greedy algorithm, so that the algorithm has good reconstruction performance no matter in a noisy environment or a non-noisy environment. But this results in a large amount of computation since atoms with different scales, frequencies, times and phases need to be screened one by one in each step to find the best atom that matches the signal structure.

Disclosure of Invention

In order to solve the above problem, the present disclosure provides a magnetic resonance image denoising method, including the following steps:

s100: inputting and initializing an image y to be reconstructed;

s200: generating an initial population W according to the upper limit and the lower limit of the parameters of atoms of the image to be reconstructed;

s300: initializing iteration times, accelerating to search L optimal atoms by using a self-adaptive genetic algorithm, respectively matching residual signals by using the L optimal atoms, updating a candidate support set and the residual signals of the atoms to obtain a candidate support set of the iteration, and repeatedly executing S300 if the iteration times are less than the maximum iteration times;

s400: carrying out secondary screening on the obtained candidate support set;

s500: when the iteration requirement is met, finding out the candidate support set with the minimum corresponding residual error, restoring the reconstructed signal, and outputting a reconstruction result graph

The genetic algorithm is simple in principle, is not limited by the problem field, is only related to the optimized objective function, and has the capability of quickly searching the problem, so that the genetic algorithm is introduced into the iterative optimization process of the multi-path matching pursuit algorithm for optimization, and the decomposition speed of the multi-path matching pursuit algorithm is increased. In order to avoid the defect of early maturity of the genetic algorithm, the method utilizes the self-adaptive genetic algorithm to automatically adjust the cross probability and the variation probability in the evolution process, overcomes the defect that the algorithm falls into local optimum, and improves the efficiency and the reconstruction precision of the algorithm.

Drawings

Fig. 1 is a schematic flow chart of a magnetic resonance image denoising method provided in an embodiment of the present disclosure;

FIGS. 2(a) to 2(f) are standard image denoising results in one embodiment of the present disclosure;

FIG. 3 is a PSNR value and runtime of a Lena graph and a Barbara graph in one embodiment of the present disclosure;

fig. 4(a) to 4(f) are simulated brain MR image denoising results in another embodiment of the present disclosure;

FIG. 5 is a PSNR value and runtime in another embodiment of the present disclosure;

6(a) to 6(e) are clinical breast MR images in one embodiment of the present disclosure;

fig. 7(a) to 7(e) are clinical brain MR images in one embodiment of the present disclosure;

fig. 8 is a breast MR image 1, a breast MR image 2, and a run time of a brain MR image in an embodiment of the present disclosure.

Detailed Description

Referring to fig. 1, in one embodiment, a magnetic resonance image denoising method is disclosed, which includes the following steps:

s100: inputting and initializing an image y to be reconstructed;

s400: carrying out secondary screening on the obtained candidate support set;

In another embodiment, step S100 further comprises:

the initialization is as follows: the number of iterations k is 0 and the initial residual r⁰Y, initial set of candidate support sets S⁰And { phi }, the maximum number of iterations Ntem.

In another embodiment, the upper and lower limits of the atomic parameter are specifically set as follows:

atom g_γThe lower limit of the parameter (x, y, s, θ) is bi ═ 0, 0, 0, 0]With the upper limit being bs ═ hx, hy, 2^j，π]Hx, hy respectively represent the height and width of the image, j represents the index of the scale factor, and j is more than or equal to 1 and less than or equal to 4; where (x, y) represents the position of the atom, s represents a scale parameter of the atom, and θ represents the deflection angle of the atom.

In another embodiment, step S300 further comprises:

s301: when k < Ntem, the number of iterations: k is k +1, candidate index: u is 0, candidate support set: s^k＝φ；

S302: when 1 is more than or equal to i and less than or equal to | S^k-1L, accelerated search of the L optimal atoms using adaptive genetic algorithm

And using the L optimal atoms by calculating proj<r_i ^k-1，g_r>Matching of residual signals is performed, wherein | S^k-1I represents the number of candidate support set sets for the (k-1) th iteration; r is_i ^k-1Represents the ith residual of the (k-1) th iteration; proj represents a variable of the calculated inner product value, i represents a loop variable;

s303: when j is more than or equal to 1 and less than or equal to L, a temporary support set is established,

wherein S is_i ^k-1An ith support set representing a k-1 th iteration candidate support set; g_rjJ-th atom representing the optimal L atoms selected for the current candidate index;

s304: checking if there is a common support set, if

The candidate index u is updated to u +1 and the current candidate support set is updated

Updating a set of candidate support sets

Then updating the reconstructed signal

Re-updating residual signals

Otherwise, returning to step S303 until j equals L, ending step S303, returning to step S302, and forming the k-th candidate support set

In another embodiment, step S400 further comprises

When the number of iteration layers reaches a certain value, namely k is more than or equal to gamma, the candidate support set S is subjected to the set matching with the residual error signal according to the size of the corresponding residual error signal^kThe candidate support set in (1) is secondarily screened to remove the residual error with the maximum value

A candidate support set, where Γ represents a constant, specifying a particular number of iteration layers; u denotes the number of candidate index sets of the current iteration layer, and α denotes a constant used to calculate the number of pruned index sets.

The method utilizes the characteristics of excellent performance and small calculated amount of the genetic algorithm to solve the multi-parameter problem, adopts the self-adaptive genetic algorithm to optimize and realize the multi-path matching tracking algorithm, overcomes the defect that the traditional genetic algorithm is easy to fall into local optimization, searches the parameters of the characteristic structure of the optimal matching image, reduces the calculated amount to a great extent, and improves the decomposition speed and the reconstruction precision of the multi-path matching tracking algorithm.

In another embodiment, a two-dimensional Gabor function is used as a generating function of the atom dictionary, the two-dimensional Gabor function is obtained by multiplying a two-dimensional gaussian window function w (x, y) by a two-dimensional complex wave s (x, y), and g (x, y) ═ s (x, y) w (x, y), where the two-dimensional gaussian window function w (x, y) and the two-dimensional complex wave s (x, y) are respectively expressed as follows:

x′＝(x-x₀)cosθ+(y-y₀)sinθ，y′＝(y-y₀)cosθ-(x-x₀)sinθ，s(x，y)＝exp(i(2π(u₀x+v₀y)) to obtain a two-dimensional Gabor kernel function, and obtaining a Gabor atomic cluster through translation and modulation, wherein the two-dimensional Gabor kernel function is expressed as follows:

wherein: x, y, x ', y' are position coordinates,

representing the ratio of the amplitudes of the Gaussian kernels, (u)₀，v₀) Representing frequency domain coordinates; four-dimensional parameter set gamma (x)₀，y₀S, θ) represents an atom, (x)₀，y₀) The position of a two-dimensional Gabor atom is represented, S represents a scale parameter of the atom, θ represents a deflection angle of the atom, and the two-dimensional complex wave S (x, y) is in a complex form.

The Gabor function proved to be similar to the receptive field of simple cells in the human visual system with good local characterization properties and directional selectivity. The method is sensitive to image edge information and has good performance in the local space domain and the frequency domain of the extracted target. Based on the advantage, the atom library generated by adopting the Gabor function as the basis function can well represent the spatial frequency (scale), the spatial position and the direction selectivity of the image.

The method takes a 4-dimensional parameter set gamma of a Gabor atom as a parameter to be optimized in the AGA (chromosome in a genetic algorithm), and takes an inner product of a residual signal and the atom as a fitness function of the genetic algorithm. The encoding mode of the parameters is real number encoding, so that no binary conversion is needed in the atom optimizing process, the algorithm is easy to understand, and the time is saved.

In another embodiment, the method calculates the probability of crossover and mutation according to the following equation. Wherein p is_{c_max}＝0.9，p_{c_min}＝0.6，p_{m_max}＝0.1，p_{m_min}The reason for this is mainly to guarantee the group at 0.01The cross probability and the variation probability of the individual with the maximum fitness value in the body are not 0, so that the algorithm jumps out of the local optimal solution.

Wherein f is_maxRepresenting the maximum fitness value in the population,

representing the mean fitness value of each generation population, f' representing the greater fitness value of the two individuals undergoing crossover, f representing the fitness value of the individual undergoing mutation, p_{c_max}Denotes the maximum crossover probability, p_{c_min}Denotes the minimum crossover probability, p_{m_max}Denotes the maximum mutation probability, p_{m_min}Representing the minimum mutation probability.

The following examples are illustrated in conjunction with fig. 2(a) to 8.

In the process of acquiring and imaging the MR image, the real part and the imaginary part of a signal are generally interfered by Gaussian additive noise with the same variance and zero mean, and the noise carried by the real part and the imaginary part is independent of each other. The obtained spatial data is subjected to a modulus operation after an inverse Fourier transform to obtain an MR image, at this time, noise in the image is no longer simple white gaussian noise, and is generally considered to be distributed in a rice (Rician) characteristic, and a probability density function is as follows:

wherein, I₀Represents a modified zero-order Bessel function of the first kind, where A represents the peak of the signal amplitude and the Rice distribution transforms into a Rayleigh distribution when A → 0. The Rice distribution can be understood as the majority of the main signal and the obedient Rayleigh distributionThe sum of the path signal components. z is the envelope of the sine (cosine) signal plus the narrow-band gaussian random signal, and o ^2 is the power of the multipath signal component.

In one embodiment, a 256 x 256 standard test chart from Lena and Barbara is used, as shown in fig. 2(a), to which is added a rice noise with a signal to noise ratio of 20, as shown in fig. 2 (b). Fig. 2(c) -fig. 2(f) show the processing results of the median filtering, MP algorithm, MMP algorithm, and the present method on fig. 2(b), respectively. Wherein the template size of the median filter is 3 x 3, and the MP algorithm, MMP algorithm, and the method all utilize 2500 atoms for reconstruction.

The denoising results were evaluated from the subjective visual and objective evaluation indices peak signal-to-noise ratio (PSNR), as shown in fig. 3. Wherein the peak signal-to-noise ratio is calculated as follows:

in the formula (I), the compound is shown in the specification,

and x_pThe restored pixel values and the pixel values of the original image, the size of the image M × N, respectively. The larger the value of the peak signal-to-noise ratio. The better the representative image recovery.

In summary, the following results can be obtained: compared with other three algorithms, the traditional median filtering denoising method has the defects of weak noise suppression, smooth result graph and edge blurring; the MP algorithm, the MMP algorithm and the algorithm in the method have better suppression on noise, the MMP algorithm has fewer noise points in an image compared with a processing result of the MP algorithm, and the method is visually different from the MMP algorithm. Compared with the traditional median filtering algorithm, the method has longer running time, but shorter than MP algorithm and MMP algorithm. The reason for this is that the MP and MMP algorithms achieve the purpose of denoising based on sparse reconstruction of the atom dictionary, and an inner product needs to be calculated for each atom in the atom dictionary library and the residual signal to obtain the best atom in the iterative process, so the time for median filtering is long; the MMP algorithm needs to select a plurality of candidate atoms in each iteration to improve the reconstruction precision, so the time is longer than that of the MP algorithm; but the method introduces the self-adaptive genetic algorithm to accelerate the decomposition and reconstruction speed of the MMP algorithm, so the time is short. Compared with the PSNR value of the median filtering and MP algorithm, the method has the advantage that the PSNR value is improved.

In another embodiment, verification is performed using the image in the database of images of simulated brain images from Mcgcill university MR. The brain MR image is a T1-weighted image, the noise level is 7%, the size is 217 × 181, and the original image and the noise map are shown in fig. 4(a) and (b), respectively. Fig. 4(c-f) are the results of median filtering, MP algorithm, MMP algorithm, and the present method, respectively, on fig. 4 (b). Wherein MP, MMP and the process all use 2000 atoms for the reconstruction.

As shown in fig. 5, the PSNR value of the median filtered image is 30.6506dB, and the PSNR values of the MP algorithm, the MMP algorithm, and the method are 30.7408dB, 30.8693dB, and 30.9597dB, respectively, which shows that compared with other three comparison algorithms, the PSNR of the method is improved by 0.3091dB, 0.2189dB, and 0.1618dB, respectively. However, compared with the MP algorithm, the MMP algorithm improves accuracy and increases the running time of the algorithm correspondingly, mainly because the MMP algorithm selects a plurality of atoms to match residual signals in an iterative process, the method optimizes the MMP algorithm by using the genetic algorithm, accelerates the algorithm decomposition and reconstruction speed, and improves the PSNR value compared with the algorithm without the optimization.

In another embodiment, two of the three clinical MR images are breast MR images, as shown in the first and third rows of fig. 6 (a). The two sets of image data were derived from a standard 1.5T dual milk coil magnetic resonance scanner, siemens, germany, as a three-dimensional T1 weighted gradient echo sequence diagram after contrast enhancement. When scanning and imaging, the patient adopts a prone position, the mammary glands on the two sides naturally hang in the cavity of the breast coil, and the parameters are as follows: TR 5.6ms, TE 2.76ms, layer spacing 0.3mm, layer thickness 1.2mm, FOV 34 × 34cm, image size 512 × 512, and scanning time 60s each; the other set is brain MR images, as shown in fig. 7 (a). The image data is from the signa HDxT 3.0T magnetic resonance imaging system of GE company and is a T1 imaging sequence chart. The human head is in the front during scanning, is the supine position, and each parameter is: TR 3042.13ms, TE 11.77ms, 6mm interlayer spacing, 5mm thickness, DFOV 220mm, and image size 512 × 512.

Fig. 6(e) and fig. 7(e) show the results of the reconstruction of two groups of data by the method, and the results of the median filtering, MP algorithm and MMP algorithm are shown for comparison, as shown in fig. 6 and fig. 7(b) - (d). For a more intuitive comparison algorithm, the region with the lesion site (tumor) in the two clinical breast MR images in fig. 6(a) was truncated (second and fourth lines in fig. 6 (a)) and compared with the comparison algorithm. Meanwhile, the run-time comparison of the algorithms is given in fig. 8.

The data according to fig. 6(a) to 6(e), 7(a) to 7(e) and 8 are comprehensively analyzed: compared with other three algorithms, the traditional median filtering method has no difference in the suppression of the noise of the whole image, but compared with a local processing result, the processed result image has the defects of fuzzy edges, incomplete detail retention and the like, and particularly in the local processing result of a second mammary gland, the smoothness degree of the image edge is large, so that certain influence is generated on the extraction of subsequent image characteristics and the like; the MP algorithm, the MMP algorithm and the method are all used for noise suppression based on sparse representation essentially, and from the processing results of fig. 6(a) to 6(e), the left breasts (with focus parts) of the two mammary gland images are seen, and the method and the MMP algorithm have better quality for gland reconstruction than the MP algorithm. In the aspect of time consumption, as can be seen from the data in fig. 8, the time of the traditional median filtering is shortest, and the MMP algorithm has a longer running time than the MP algorithm, mainly because the MMP algorithm is an improvement on one candidate support set of the MP algorithm, so that a plurality of candidate support sets exist after each iteration is completed, but the method utilizes AGA to accelerate MMP decomposition and reconstruction. Therefore, the method can well retain the effective information of the original image, well inhibit noise, and relatively not increase time consumption, and can play a certain auxiliary role in subsequent clinical treatment.

Although the embodiments of the present invention have been described above with reference to the accompanying drawings, the present invention is not limited to the above-described embodiments and application fields, and the above-described embodiments are illustrative, instructive, and not restrictive. Those skilled in the art, having the benefit of this disclosure, may effect numerous modifications thereto without departing from the scope of the invention as defined by the appended claims.

Claims

1. A magnetic resonance image denoising method, comprising the steps:

S100: Input the image to be reconstructed y and initialize it;

S200: Generate an initial population W according to the upper and lower limits of the parameters of the atoms of the image to be reconstructed;

S300: Initialize the number of iterations, use the adaptive genetic algorithm to speed up the search for L optimal atoms, and use the L optimal atoms to match the residual signals respectively, update the atomic candidate support set set and the residual signal, and obtain this The candidate support set set of the next iteration, if the number of iterations is less than the maximum number of iterations, repeat S300;

S400: Perform secondary screening on the obtained set of candidate support sets;

S500: When the iteration requirement is met, find the candidate support set set with the smallest corresponding residual, restore the reconstructed signal, and output the reconstruction result graph

Step S100 further includes:

The initialization is as follows: the number of iterations k=0, the initial residual r ⁰ =y, the initial candidate support set S ⁰ ={φ}, the maximum number of iterations Ntem;

The upper and lower limits of the parameters of the atoms are specifically set as follows:

The lower limit of the parameters of atom g _γ = (x, y, s, θ) is bi = [0, 0, 0, 0], and the upper limit is bs = [hx, hy, 2 ^j , π], hx, hy represent Image height and width, j represents the index of the scale factor, 1≤j≤4; where (x, y) represents the position of the atom, s represents the scale parameter of the atom, and θ represents the deflection angle of the atom;

Step S300 further includes:

S301: when k<Ntem, the number of iterations: k=k+1, the candidate index: u=0, the candidate support set set: ^Sk =φ;

S302: When 1≤i≤|S ^k-1 |, use adaptive genetic algorithm to speed up the search for L optimal atoms

And use these L optimal atoms to match the residual signal by calculating proj=< _r _i ^k-1 , gr >, where |S ^k-1 | represents the candidate support set for the k-1th iteration The number of sets; ri _k ^-1 represents the ith residual of the k-1th iteration; proj represents the variable of the calculated inner product value, and i represents the loop variable;

S303: When 1≤j≤L, establish a temporary support set,

In Ju, S _i ^k-1 represents the i-th support set of the k-1 th iteration candidate support set set; g _rj represents the j-th optimal L atom selected for the current candidate index;

S304: Check whether there is a common support set, if

Then update the candidate index u=u+1, update the current candidate support set

Update the set of candidate support sets

Then update the reconstructed signal

Update the residual signal again

Otherwise, go back to step S303, until j=L, end step S303, go back to step S302, and form the kth candidate support set set

2. A magnetic resonance image denoising method according to claim 1, wherein step S400 further comprises:

When the number of iterative layers reaches a certain value, that is, when k≥Γ, according to the size of the corresponding residual signal, the candidate support sets in the candidate support set set S ^k are screened twice, and the one with the largest residual error is removed.

A candidate support set, where Γ represents a constant used to specify a specific number of iteration layers; u represents the number of candidate index sets for the current iteration layer, and α represents a constant used to calculate the number of index sets to be pruned.

3. a kind of magnetic resonance image denoising method according to claim 1, the generation function of the atom in step S200 is two-dimensional Gabor function, and two-dimensional Gabor function g (x, y) is by two-dimensional Gaussian window function w (x, y) is obtained by multiplying the two-dimensional complex wave s(x, y), g(x, y) = s(x, y) w(x, y), where the two-dimensional Gaussian window function w(x, y) and the two-dimensional complex wave s(x, y) are respectively expressed as follows:

x'=(xx ₀ )cosθ+(yy ₀ )sinθ, y'=(yy ₀ )cosθ-(xx ₀ )sinθ, s(x,y)=exp(i(2π(u ₀ x+v ₀ y ))), the two-dimensional Gabor kernel function is obtained by arranging the above formula, and the Gabor atomic cluster is obtained by translation and modulation, which is expressed as follows:

Where: x, y, x', y' are the position coordinates,

represents the ratio of the Gaussian kernel amplitude, (u ₀ , v ₀ ) represents the frequency domain coordinates; the four-dimensional parameter group γ(x ₀ , y ₀ , s, θ) represents an atom, and (x ₀ , y ₀ ) represents the two-dimensional Gabor atom The position of , s represents the scale parameter of the atom, θ represents the deflection angle of the atom, and the two-dimensional complex wave s(x, y) is a complex number.

4. a kind of magnetic resonance image denoising method according to claim 1, in order to speed up the speed of searching L optimal atoms in step S300, in adaptive genetic algorithm, determine crossover probability and mutation probability according to the following formula,

Among them, _fmax represents the maximum fitness value in the group,

Represents the average fitness value of each generation of the population, f' represents the larger fitness value of the two individuals performing crossover, f represents the fitness value of the individual performing mutation, p _{c_max} represents the maximum crossover probability, and p _{c_min} represents the minimum crossover Probability, p _{m_max} represents the maximum mutation probability, p _{m_min} represents the minimum mutation probability; wherein, set p _{c_max} = 0.9, p _{c_min} = 0.6, p _{m_max} = 0.1, p _{m_min} = 0.01, to ensure the crossover probability of the individual with the largest fitness value in the population and the mutation probability is not 0, which makes the algorithm jump out of the local optimal solution.