CN104751423B

CN104751423B - A kind of Medical Image Denoising method based on resonance subspace analysis

Info

Publication number: CN104751423B
Application number: CN201510114525.2A
Authority: CN
Inventors: 梁军利; 柯婷; 于国阳; 贾薇; 叶欣; 范文; 李敏
Original assignee: Xian University of Technology
Current assignee: Xian University of Technology
Priority date: 2015-03-16
Filing date: 2015-03-16
Publication date: 2017-09-29
Anticipated expiration: 2035-03-16
Also published as: CN104751423A

Abstract

The invention discloses a kind of Medical Image Denoising method based on resonance subspace analysis, in units of similar block group, initial filter model is built first, and initial filter is carried out to image, then similar image block is utilized in the resonance characteristics of transformation space, add and protect side driving and the holding of local manifolds structure setting, improve Filtering Model, further image is repeatedly filtered, filtering is based on last time filtered image each time, so that the image block after denoising keeps image edge information to greatest extent, so as to farthest recover image information；The denoising method of the present invention can effectively reduce the noise of medical image so that the medical image after noise reduction has the when superior visual effect of good noise, and the image after obtained denoising can provide good basis for follow-up image procossing.

Description

A kind of Medical Image Denoising method based on resonance subspace analysis

Technical field

The invention belongs to image processing field, it is related to a kind of Medical Image Denoising method based on resonance subspace analysis.

Background technology

Medical Imaging Technology is important part in present medical science, and as one with fastest developing speed in medical technology Technology.Medical Imaging Technology mainly includes medical imaging Display Technique, medical image analysis treatment technology and medical image pressure Three Main ways of contracting transmission technology.Its main function is：Gather the information of patient body lesion portion and be stored as corresponding Image, passes through the state of an illness clear, in detail for learning patient to these image informations are for further analysis, diagnosis comes.So as to more Good carries out further treatment to patient.The image information of reservation is also used as the reference diagnosed in the future.Modern medicine shadow As technology also makes it possible tele-medicine, the communication between doctor and patient is greatly facilitated.

Because the principle of imaging is different with equipment, medical imaging technology has a variety of imaging patterns.For in terms of big, The anatomy imaging pattern of physiologic form and the functional imaging pattern of description bodily fuctions or metabolism, main bag can be allocated as describing Include：X-ray imaging, Computed tomography (CT), magnetic resonance imaging (MRI), ultrasonic imaging (US), positron emission are broken Layer scanning imagery (PET), electroencephalogram (EEG).

Using digital image processing techniques and computer technology to rely on, the analysis and processing of medical image are medical image skills A particularly important link in art, it be doctor obtain patient's authentic data information important guarantee, be also doctor carry out into The important evidence of one step treatment.The analyzing and processing of medical image, mainly includes：The pretreatment of medical image, feature extraction, medical science Image registration, medical image enhancement, medical image segmentation, Medical Images Classification and Medical Image Compression, storage are with communicating Deng.

With the development of technology, medical imaging technology is also being improved constantly, and medical image quality is also improved therewith.However, Medical image in imaging process, due to medical image acquisition, transmission and obtain during always inevitably by The interference of various noise sources, these noises may be interfered to the diagnosis of doctor.Thus, corresponding pre- place is done to medical image The noise included in reason, reduction image, improves signal to noise ratio, is extremely important.

Medical image filtering is a basic and crucial technology in Medical Image Processing, is always at medical image One problem in reason field.People are according to the characteristics of practical medical image, the regularity of distribution of the statistical property of noise and frequency spectrum, carry Many denoising methods are gone out.Two major classes can be mainly divided into：One class is the filtering method based on spatial domain, i.e., in pending image Pointwise movement spatial filter (that is, spatial domain template), wave filter responded in the point by the filter coefficient of predefined with The relation of the respective pixel values in Filtering Template is inswept region is calculated, and it represents algorithm mean filter, non-local mean filtering (NLM) etc.；It is another kind of to be referred to as frequency domain filtering, pending image is exactly transformed into frequency domain, by the processing to frequency coefficient, then Inverse transformation obtains filtering image, and it, which represents algorithm, wavelet transformation, singular value decomposition etc..

But above-mentioned filtering method is all had some limitations and unreliability, and image has been obscured to a certain extent What is more so that image produces distortion, the picture quality after denoising is low, is that the processing of follow-up medical image is caused to a certain degree Difficulty.Therefore, advanced, the effective denoising method of medical image is very important, and this is related at follow-up medical image Reason or the progress of medical diagnosis.

The content of the invention

It is an object of the invention to provide a kind of Medical Image Denoising method based on resonance subspace analysis, solve existing The low technical problem of medical image quality in technology after denoising.

The technical solution adopted in the present invention is that a kind of Medical Image Denoising method based on resonance subspace analysis has Body is implemented according to following steps：

Step 1：The medical image that a width contains random noise is inputted, this noisy medical image is designated as Y, Y ∈ R^M×N, Y=X+ E, X are original medical image, and E is the white Gaussian noise of addition, and M × N represents the size of image；

Step 2：Y is designated as to noisy medical image and carries out initial filter, detailed process is as follows：

2.1, noisy medical image Y is divided into image block of several sizes as n × n in the way of sliding block step-length is 1；

2.2, noisy medical image Y noise criteria difference σ is estimated with the robustness mediant estimation method based on small echo threshold₁；

2.3, it is target image block that each image block is selected successively, and its similar block is screened using threshold method, builds similar block group；

2.4, calculate the weighted value of the image block in similar block group；

2.5, in units of similar block group, initial filter is carried out to the target image block in each similar block group respectively, specifically Method is as follows：

2.5.1, the initial filter model that the characteristic for noisy medical image Y is built is as follows：

Constraints

Wherein, η, γ are preset constant, K₁Represent the image number of blocks included in similar block group, P_kRepresent in similar block group K-th of image block,H is denoising parameter, U₁={ u₁,u₂,...u_d},U₁∈R^n×d, V₁={ v₁, v₂,...v_d},V₁∈R^n×d, U₁,V₁For left and right basic matrix, it is made up of orthogonal base vectors, f₁Estimate for the filtering of image block,Generation The weighted value of each image block in table similar block group,Square of F norms is represented, | | | |_TVRepresent TV norms；

2.5.2, with alternating iteration multiplier method, the U that initial filter model corresponds to each target image block is solved₁,V₁；

2.5.3, the U obtained according to solution₁,V₁, obtain estimating corresponding to the filtering of each target image block：

2.5.4, the filter result of each target image block is added to the one and noisy medical image Y blank square with size On the corresponding position of battle array, the number average value obtained after each filtering for taking each pixel as this pixel last estimation Value, obtains the initial filter result to noisy medical image；

Step 3：Using the initial filter result of noisy medical image as initialization, further filtering, specific method is as follows：

3.1, according to the method for step 2.1~step 2.4 successively by last time filtered noisy medical image piecemeal, estimation Its noise criteria is poor, builds similar block group, calculates the weighted value of each image block in similar block group；

3.2, in units of similar block group, the target image block in each similar block group is filtered again respectively, specifically Method is as follows：

3.2.1, the image denoising model based on resonance subspace analysis is built：

Constraints

Wherein, K₂Represent the image number of blocks included in similar block group, M_kRepresent k-th of image block, M in similar block group_k∈ R^n×n, U₂={ u '₁,u′₂,...u′_d},U₂∈R^n×d, V₂={ v '₁,v′₂,...v′_d},V₂∈R^n×d, U₂,V₂For left and right basic matrix, It is made up of orthogonal base vectors,H is denoising parameter, f₂Estimate for the filtering of image block,Represent The matrix of the composition reciprocal of pixel filter times in image block,J₀Represent with last filtered figure The image block divided as based on, that is, the pending target image block of this filtering,Represent in similar block group The weighted value of each image block；

3.2.2, with alternating iteration multiplier method, the denoising model in solution procedure 3.2.1 corresponds to the U of each image block₂、 V₂；

3.2.3, the U obtained according to solution₂、V₂, obtain corresponding to the filtering estimation for respectively obtaining image block：

3.2.4, the filter result of each target image block is added to the one and noisy medical image Y blank square with size On the corresponding position of battle array, the number average value obtained after each filtering for taking each pixel as this pixel last estimation Value, obtains the filter result again to noisy medical image；

Repeat repeatedly to filter image using the filtering method of step 3, after filtering each time is with last time filtering Based on obtained image, each image block or so basic matrix is calculated to being filtered to image, until left and right basic matrix is to convergence, Obtain the medical image figure after final denoising.

The features of the present invention is also resided in,

The calculation formula of noise criteria difference is as follows in step 2：

σ₁=median (W_HH)/0.6745；

Wherein, median (W_HH) represent that image decomposes through 2-d wavelet after, the wavelet systems of obtained high pass-high pass subspace Number.

The calculation formula of the weighted value of each image block in step 2 in similar block group is as follows：

α represents time number formulary in above formula.

The calculation formula of the weighted value of each image block in step 3 in similar block group is as follows：

α represents time number formulary in above formula.

The method that similar block group is built in step 2 is as follows：

A, first selected digital image block are used as target image block P₀, P₀∈R^n×n；

B, then sets distance threshold τ_d=3 σ₁ ²n², successively to be sized in target image block centered on each pixel For x × x search window N (s), the other image blocks and target image block P in the range of each search window are calculated₀Euclidean distance, calculate Formula is as follows：

Wherein, Y (t) is the gray value of other image blocks, and Y (s) is the gray value of target image block；

C, finally judges and target image block P₀Similar image block, deterministic process is as follows：

If d (s, t)≤τ_d, then the image block in search window and target image block P₀It is similar, P is designated as successively₁、P₂……、If d (s, t) ＞ τ_d, then image block and target image block P₀It is dissimilar；

D, similar block group is built into by the similar image block selectedP_k∈R^n×n, k=0,1,2 ..., K₁；

U is solved in step 2.5.2₁,V₁Method it is as follows：

2.5.2.1, Lagrange multiplier λ is introduced₁, relaxation object function bound term：

Wherein, λ₁It is Lagrange multiplier, its dimension and f₁Unanimously, μ is preset constant；

2.5.2.2, each variable of initial filter model is initialized, U is initialized using singular value decomposition SVD₁、V₁ForV₁ ¹, λ₁The null matrix being initialized as with image block formed objectsf₁Target image block is initialized as, is designated as

2.5.2.3, solution by iterative method U₁、V₁, iterations is designated as g；

(1) f is solved₁ ^g+1：Solve f₁Method it is as follows：

Gradient descent method updates f₁Formula be：, z represent gradient descent method update f₁Iteration time Number, dt represents to spread step-length, f₁₍₁₎=f₁ ^g；

Simultaneous above formula：WillV₁=V₁ ^gSubstitution formula produces f₁ ^g+1；

(2) alternating iteration Optimization SolutionV₁ ^g+1：

1. U is initialized₁,V₁For orthogonal matrix, it is designated as

2. solveM is alternating iteration Optimization SolutionIterations：

Fixed V₁, solve U₁Method it is as follows：

To formulaEigenvalues Decomposition, By the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value, that is, constitute U before taking₁；

Willf₁=f₁ ^g+1,Substitute into, obtain

3. solveFixed U₁, solve V₁Method it is as follows：

To formulaCarry out characteristic value Decompose, by the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value, as V before taking₁；

Willf₁=f₁ ^g+1,Substitute into, obtain

Wherein,L₁ =D₁-S₁,I is unit battle array,

4., withFor next iteration process U₁,V₁Initial value, repeat step 2., 3., until convergence, ObtainV₁ ^g+1；

(3) obtained using above-mentioned stepsV₁ ^g+1、f₁ ^g+1CalculateCalculation formula is as follows：

(4) withV₁ ^g+1、f₁ ^g+1、For next iteration process U₁、V₁、f₁、λ₁Initial value, repeat step (1) ~(3), until the convergence of each variable, obtains final U₁,V₁。

U is solved in step 3.2.2₂、V₂Method it is as follows：

3.2.2.1, Lagrange multiplier λ is introduced₂, relaxation object function bound term：

Wherein, μ is default constant；

3.2.2.2, each variable of denoising model in initialization step 3.2, U is initialized using singular value decomposition SVD₂,V₂ForBy λ₂It is initialized as the null matrix with image block formed objectsf₂It is initialized as the similar of last filtered image Target image block in block group, is designated as

3.2.2.3, solution by iterative method U₂、V₂, iterations is designated as c：

3.2.2.3.1, solveSolve f₂Method it is as follows：

Gradient descent method updates f₂Formula be：L represents that gradient descent method updates f₂Iteration Number of times, dt represents diffusion step-length,

Simultaneous above formula：WillSubstitution formula is produced

3.2.2.3.2, solve

(1) U is initialized₂,V₂For orthogonal matrix, it is designated asξ be initialized as with tile size identical null matrix, It is designated as ξ¹, alternately solve t, V₂、U₂, ξ iterations be designated as d；

(2) t is sought^d：

U is fixed using least square method₂、V₂Solve t method be：

By expression formulaCalculated using least square method and obtain T, then this is vectorial The matrix of corresponding form is converted to, t is produced；

Wherein, F, T, R, b₂Respectively matrixt,r,b₁Vectorization form, B_bFor diagonal matrix, element is vector b on its diagonal₂The element of middle correspondence position, b₁For target image block The matrix of correspondence position, b in B matrixes₁∈R^n×n, B matrixes are the corresponding filter times composition of entire image each pixel Matrix；

Willξ=ξ^dSubstitute into, produce t^d；

(3) solveSolve first

A, initializes U₂,V₂For orthogonal matrix, it is designated as in step 3.2.2.3.2 iterative process

B, is solvedP is alternative method solution in step 3.2.2.3.2Iteration Number of times：

1. solveFixed V₂Seek U₂Method it is as follows：

To formulaCarry out feature Value is decomposed, by the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value before taking, and constitutes U₂；

By t=t^d,ξ=ξ^dSubstitute into, obtain

2. solveFixed U₂Seek V₂Method it is as follows：

To formulaCarry out special Value indicative is decomposed, by the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value before taking, and constitutes V₂；

By t=t^d,ξ=ξ^dSubstitute into, obtain

Wherein,ξ is Lagrange multiplier, and κ is preset constant, L₂=D₂-S₂, I is unit battle array, i=1,2 ..., K₂：

C, withFor the U of next iteration process₂、V₂Initial value, repeat step 1., 2. until receive Hold back, obtain

(4) ξ is solved^d+1：

(5) withξ^d+1,For the U of next iteration process₂、V₂、ξ、λ₂、f₂Initial value, Repeat step (2)~(4), until convergence, is obtained

3.2.2.3.3, obtained using above-mentioned stepsCalculateCalculation formula is as follows：

3.2.2.3.4, calculate what is obtained with above-mentioned stepsDuring next iteration U₂、V₂、f₂、λ₂Initial value, repeat step 3.2.2.3.1~3.2.2.3.3 until each variable restrain, obtain final U₂, V₂。

The beneficial effects of the invention are as follows the present invention, in the resonance characteristics of transformation space, is added using similar image block and protects side Driving and the holding of local manifolds structure are set, and one group of optimal basic matrix are trained, for constructing filtering image block so that go Image block after making an uproar keeps image edge information to greatest extent, so as to farthest recover image information；The present invention's goes Method for de-noising can effectively reduce the noise of medical image so that the medical image after noise reduction has good noise is when superior to regard Feel effect, the image after obtained denoising can provide good basis for follow-up image procossing.

Brief description of the drawings

Fig. 1 is the noisy medical image of noise variance σ=20；

Fig. 2 is after being filtered in the prior art using mean filter method to the noisy medical image of noise variance σ=20 Result figure；

Fig. 3 is after being filtered in the prior art using median filtering method to the noisy medical image of noise variance σ=20 Result figure；

Fig. 4 is the result figure after being filtered using the method for the present invention to the noisy medical image of noise variance σ=20；

Fig. 5 is the noisy medical image of noise variance σ=30；

Fig. 6 is after being filtered in the prior art using mean filter method to the noisy medical image of noise variance σ=30 Result figure；

Fig. 7 is after being filtered in the prior art using median filtering method to the noisy medical image of noise variance σ=30 Result figure；

Fig. 8 is the result figure after being filtered using the method for the present invention to the noisy medical image of noise variance σ=30.

Embodiment

Resonance subspace (Concurrent Subspace) is the member for having similar characteristic by some in transformation space Show the subspace of resonance characteristics.

A kind of Medical Image Denoising method based on resonance subspace analysis of the present invention, implements according to following steps：

Step 1：The medical image that a width contains random noise is inputted, this noisy medical image is designated as Y, Y ∈ R^M×N, Y=X+ E, X are original medical image, and E is the random noise of addition, and M × N represents the size of image, the present invention handle for Gauss white noise Sound；

Step 2：Initial filter is carried out to noisy medical image Y, detailed process is as follows：

2.1, by the noisy medical image Y piecemeals in step 1, by noisy medical image with sliding block in present embodiment The mode that step-length is 1 is divided into the image block that several sizes are n × n, n=8.

2.2：Estimate that noisy medical image Y noise criteria is poor, with the robustness mediant estimation method based on small echo threshold (Median Absolute Deviations, MAD) carries out noise estimation, and its mathematic(al) representation is as follows：

σ₁=median (W_HH)/0.6745 (1)；

Wherein, median (W_HH) represent after being decomposed through 2-d wavelet, the wavelet coefficient of obtained high pass-high pass subspace； The theoretical foundation of this method：Noise wavelet coefficients are concentrated mainly on high-frequency sub-band, so noise criteria difference can be by high frequency certainly The coefficient of band is approximate；

2.3, it is target image block that each image block is selected according to this, and the similar block of target image block, structure are screened using threshold method Build similar block group：

2.3.1, selected target image block P first₀, P₀∈R^n×n；

2.3.2, then distance threshold τ is set_d=3 σ₁ ²n², set successively centered on each pixel in target image block Size is 20 × 20 search window N (s), calculates other image blocks and target image block P in the range of each search window₀Euclidean away from From calculation formula is as follows：

2.3.3, finally judge and target image block P₀Similar image block, deterministic process is as follows：

2.3.4, the similar image block selected is built into similar block groupP_k∈R^n×n, k=0,1,2 ..., K₁；

In order that denoising effect is excellent, τ is chosen_dPrinciple it is as follows：The Gauss white noise added in known original medical image Sound submits to N (0, σ) distributions, we can assume that similar block and target image block come from same clean image block, due to by Random noise influences, and performance slightly has difference, then following stochastic variable：Obey χ²(n²) distribution, χ²(n²) The accumulation of stochastic variable can be byRepresent, (x a) is defined as wherein incomplete gamma functions γ For complete gamma function.As z >=3, for any x >= 3z, there is F (x, z) >=0.99, thus for image block and given noise criteria difference σ that size is n × n, we can choose Distance threshold τ_d=6 σ²n², we assert that image block has very high similitude in this threshold range.But actual emulation knot Fruit proves, may not be similar to target image block with the image block that this threshold value is chosen, some image blocks and object block in structure simultaneously Inconsistent, this threshold value is not reasonable, and threshold value is excessive to be caused to falsely drop result.In order to eliminate error, it is observed that as Y (t) and Y (s) to come from same clean image block, Y (s)-Y (t) value belongs toExperiments verify that, work as τ_d=3 σ²n² When, it can effectively choose suitable similar block；

By assuming that the similar block and target image block that find are to the addition of difference on the basis of same clean image block The noise of intensity and obtain, then these similar blocks should have resonance effects on projector space, thus can pass through resonon Spatial analysis, obtains clean image；

2.4, the weighted value of all image blocks in each similar block group is calculated, formula is as follows：

α represents time number formulary in above formula；

Constraints

Wherein, η, γ are η=10, γ=0.5 in preset constant, present embodiment.U₁={ u₁,u₂,...u_d, U₁∈ R^n×d, V₁={ v₁,v₂,...v_d, V₁∈R^n×d, U₁,V₁For left and right basic matrix, it is made up of orthogonal base vectors, Parameter h=30000, f₁Estimate for the filtering of target image block,Square of F norms is represented, | | | | TV represents TV norms, TV norms are defined as follows

u_x, u_yPartial differential is shown as in continuous domain, in discrete domain then shows as horizontal and vertical difference, digital picture Pixel can be considered discrete point, thus, it is assumed that u is a pixel on image, then its TV norm, is its horizontal difference and longitudinal direction The evolution of the quadratic sum of difference.

Every meaning and effect in initial filter model：Section 1 requires that the left and right basic matrix that training is obtained can cause weight Each filtering image block error sum in similar block group after structure is minimum, i.e., reconstructed error is minimized；Section 2 is required after filtering Target image block gradient minimisation, because the randomness of noise and scalar property cause the gradient of image to increase, thus add In this canonical constrain, the purpose of noise reduction can be played since while be that gradient is limited, can also optimize image border, Image detail information is kept, the effect of guarantor side is played；Section 3 is then to construct neighbour's figure according to LAP, it is desirable to the group moment that training is obtained It is neighbouring on transformation space that battle array meets similar image block.

2.5.2, in initial filter model each variable solution procedure：

With alternating iteration multiplier method (ADMM), to the known variables U in initial filter model (4)₁、V₁It is iterated excellent Change；

Wherein, λ₁It is Lagrange multiplier, its dimension and f₁Unanimously, μ is preset constant, is set in the embodiment For μ=1；

2.5.2.2, initial filter model (4) each variable is initialized, U is initialized using singular value decomposition SVD₁、V₁Forλ₁The null matrix being initialized as with image block formed objectsf₁Target image block is initialized as, is designated as

2.5.2.3, solution by iterative method U₁、V₁, solution by iterative method U₁、V₁Iterations be designated as g：

(1) f is solved₁ ^g+1：Solve f₁Method it is as follows：

Fixed λ₁,U₁, V₁, taken out and f from (6) formula₁Related item

Merge abbreviation to (7) formula to obtain

Wherein,

Solution f can be seen that by (8) formula₁The problem of can be solved with the method for solving of partial differential equation (PDE), i.e., should Optimized with gradient descent method, two in formula (8) can be considered two kinds of constraints of TV norm minimums, that is, ensure image block Under conditions of fidelity, TV norms are minimized, the influence that noise is caused to gradient is eliminated, recover image gradient information, also Recover the detailed information such as the Edge texture of image.

Then an energy functional is defined：

Then solve f₁The problem of be converted into the problem of functional seeks extreme value, i.e. variation (TV) problem；

The diffusion equation of TV Restoration models is as follows：

Gradient descent method updates f₁Formula be：

Wherein, z represents that gradient descent method updates f₁Iterations, z=1,2 ..., 80, dt represent spread step-length, dt =0.2, f₁₍₁₎=f₁ ^g；

WillV₁=V₁ ^gSubstitution formula (13) produces f₁ ^g+1；

(2) solveSolve U₁,V₁Method it is as follows：

Fixed f₁, λ₁, with U in taking-out type (6)₁,V₁Related item：

Arrangement above formula abbreviation is of equal value to be obtained：

Wherein,

This step needs to solve U₁,V₁Two variables, U₁,V₁Between be related, thus need use iteration alternative optimization Mode, that is, fix U₁Seek V₁, fixed V₁Seek U₁, continuous alternating iteration optimization, until convergence；

Alternating iteration Optimization Solution U₁,V₁The step of it is as follows：

A, initializes U₁,V₁For orthogonal matrix, it is designated asAlternating iteration Optimization SolutionIterations It is designated as m；

B, is solvedFixed V₁, solve U₁Method it is as follows：

The formal expansion of trace of a matrix will be pressed on the right of (15) formula equation, can obtained

Expansion abbreviation is arranged：

Wherein in formula (17) Section 2 byExpansionization Letter and obtain, specific simplification process is as follows：

OrderI is unit battle array,

Wherein, D_iiFor Laplacian (LAP) weight matrix each column weights sum, a square formation is this is defined herein as, i.e.,I is unit battle array,

By (17), formula is obtained：

By to formula： Eigenvalues Decomposition, by the descending sequence of characteristic value, take the corresponding characteristic vector of preceding d (d=4) individual characteristic value, that is, constitute U₁；

Willf₁=f₁ ^g+1,Substitute into, obtain

C, is solvedFixed U₁, solve V₁Method it is as follows：

With fixing V in step b₁, solve U₁Derivation method it is identical, can obtain

Wherein,

By to formula： Eigenvalues Decomposition, by the descending sequence of characteristic value, take the corresponding characteristic vector of preceding d (d=4) individual characteristic value, as V₁；

Willf₁=f₁ ^g+1,Substitute into, obtain

D, withFor next iteration process U₁,V₁Initial value, repeat step b, c, until convergence, obtain ArriveV₁ ^g+1；

(4) withV₁ ^g+1、f₁ ^g+1、For next iteration process U₁、V₁、f₁、λ₁Initial value, repeat step (1)~ (3), until each variable is restrained, final U is obtained₁,V₁；

2.5.3, the U obtained according to solution₁,V₁, obtain the filtering estimation of target image block：

2.5.4, the filter result of each image block is added to the one and noisy medical image Y blank matrix phase with size On the position answered, the number average value obtained after each filtering for taking each pixel is obtained as the last estimate of this pixel To the initial filter result to noisy medical image.

By the filtering of said process, the noise of medical image is effectively suppressed, and image detail texture and side Edge information etc. is also able to keep well.

Step 3：Using the initial filter result of noisy medical image as initialization, further filtering：

3.2, in units of similar block group, the target image block in each similar block group is filtered respectively, specific method It is as follows：

Constraints

Wherein, η, γ are η=10, γ=0.5, M in preset constant, present embodiment_k∈R^n×n, represent similar block K-th of image block in group, K₂Represent the image number of blocks included in similar block group, U₂={ u'₁,u'₂,...u'_d},U₂∈R^n×d, V₂={ v'₁,v'₂,...v'_d},V₂∈R^n×d, U₂,V₂For left and right basic matrix, it is made up of orthogonal base vectors,Parameter h=30000, f₂Estimate for the filtering of image block,Represent each figure in similar block group As the weighted value of block,Square of F norms is represented, | | | |_TVRepresent TV norms；

Expression formula in constraints：Every physical significance it is as follows：

The matrix that B is the corresponding filter times composition of entire image each pixel is defined, because the present invention be directed to have Overlapping image block is handled, so in addition to the pixel at four most edges of image, each pixel not only passes through one Sliding block step-length is 1 in secondary filtering process, the present invention, thus B matrix forms are as follows：

B=[b₁₁,b₁₂,...,b_1N；......；b_M1,b_M2,...,b_MN]

Remember b₁For the matrix of currently pending image block correspondence position in B matrixes, b₁∈R^n×n, the list of elements in matrix Show in this image block in the corresponding filter times of each pixel, modelIt represents the square of the composition reciprocal of filter times Battle array, then in modelPhysical significance be：The filtering estimation that this filtering process is obtainedFinal contribution to image block is only accounted forThis part of ratio；

A is defined as follows：

Wherein, J₀Represent the image block divided based on last filtered image, that is, this filtering Pending target image block, then A physical significance is effect of the filter result to this pending image block before retaining, this Filtering is carried out on the basis of being filtered in last time, is updated from there through continuous iteration, and the final whole structure for causing image reaches To optimal.

3.2.2, with alternating iteration multiplier method (ADMM), solve denoising model (23) and correspond to each target image block U₂、V₂；

Wherein, μ is to be set to μ=0.5 in default constant, present embodiment.

3.2.2.2, each variable in initialization denoising model (23), U is initialized using singular value decomposition SVD₂,V₂ForBy λ₂It is initialized as the null matrix with image block formed objectsf₂Target image block in initial similar block group, note For

3.2.2.3.1, solveSolve f₂Method it is as follows：

Taken out and f from (24) formula₂Related item

Merge abbreviation to (25) formula to obtain

Wherein,

Then an energy functional is defined：

Gradient descent method updates f₂Formula be：

L represents that gradient descent method updates f in formula₂Iterations, l=1,2 ..., 80, dt represent spread step-length, dt= 0.2,

Simultaneous above formula, willSubstitution is produced

3.2.2.3.2, solveSolve U₂,V₂Method it is as follows：

With U in taking-out type (24)₂,V₂Related item：

Wherein,

Due in formula exist withDot product problem, it is difficult to direct matrix is solved, thus is introduced back into variableSimplify solution procedure：

Formula (32) is equivalent to

Constraints

The solution of this link is carried out using ADMM methods again, similarly, a Lagrange multiplier ξ is introduced, it is first right Formula (33) does loose constraint：

In formula, κ is preset constant, κ=1；

The solution procedure of formula (34) is also the alternative optimization mode that uses, and iteration is carried out, that is, fixes t and seek U₂,V₂, fixed U₂,V₂ T is sought, iteration is optimized successively, until restraining to obtain U₂,V₂Solution；

(2) t is sought^d：Fixed U₂,V₂The method for seeking t is as follows：

With t continuous items in (34) formula of taking-up：

Wherein,

For the ease of solving, each matrix in formula (35) is pulled into vector and solved, this is a replacement of equal value：

Wherein, F, T, R, b₂Respectively matrixt,r,b₁Vectorization form, B_bFor diagonal matrix, element on its diagonal For vectorial b₂The element of middle correspondence position；Thus, the problem of solving t can be converted into the least square problem of vector T, its most young waiter in a wineshop or an inn The form of expression for multiplying solution is：

T is obtained by above-mentioned expression formula, then this vector is converted into the matrix of corresponding form, t is produced；

Willξ=ξ^dSubstitute into, produce t^d；

(3) solveSolve first

Solve U₂,V₂Method it is as follows：With { U in (34) formula of taking-up₂,V₂Continuous item

Wherein,

1. solveFixed V₂Seek U₂Method it is as follows：The formal expansion of trace of a matrix will be pressed on the right of (38) formula equation, It can obtain：

OrderI be unit battle array, i=1,2 ..., K₂：, I is unit battle array,Obtained using initial filter identical simplifying method abbreviation (39)：

Thus, by formula：It is special Value indicative is decomposed, and by the descending sequence of characteristic value, takes the corresponding characteristic vector of preceding d (d=4) individual characteristic value, you can constitute U₂；

By t=t^d,ξ=ξ^dSubstitute into, obtain

2. solveFixed U₂Seek V₂Method it is as follows：

Wherein,

To formula： Eigenvalues Decomposition is carried out, by the descending sequence of characteristic value, the corresponding characteristic vector of preceding d (d=4) individual characteristic value is taken, you can structure Into V₂；

By t=t^d,ξ=ξ^dSubstitute into, obtain

(4) ξ is solved^d+1：

3.2.2.3.4, calculate what is obtained with above-mentioned stepsDuring next iteration U₂、V₂、f₂、λ₂Initial value, repeat step 3.2.2.3.1~3.2.2.3.3 until each variable restrain, obtain final U₂, V₂；

3.2.4. the filter result of each image block is added to the one and noisy medical image Y blank matrix phase with size On the position answered, the number average value obtained after each filtering for taking each pixel is obtained as the last estimate of this pixel To the filter result again to noisy medical image.

In whole image filtering, the structure on model is extremely important, and invention introduces non local Thought, makes full use of the characteristic that there is similar image block in image, constructs weight, is weighted polymerization, and be currently pending Image block adds available information；Meanwhile, in order to avoid there is image blurring, the result of edge-smoothing, the present invention absorbs partially micro- Divide the guarantor side characteristic of equation, introduce full variation (TV) item, the guarantor side for model drives；In order to determine the holding of local message Property, and subspace stability, introduce manifold structure, it is ensured that similar image block is close transformation space, add calculation The stability of method.

It is the related data using denoising method of the present invention to noisy medical image filtering result below：

The noisy medical image filtering result of the noise criteria of table 1 difference σ=20

	Abdominal cavity (CT)	Thoracic cavity (CT)	Chest (CT)
				PSNR	32.8643	31.2350	33.7564
SSIM	0.9133	0.8759	0.9104

Note：Noisy image PSNR=22.1129, SSIM=0.3689

The noisy medical image filtering result of the noise criteria of table 2 difference σ=30

	Abdominal cavity (CT)	Thoracic cavity (CT)	Chest (CT)
				PSNR	30.6529	29.1996	31.6594
SSIM	0.8540	0.8174	0.8415

Note：Noisy image PSNR=18.5924, SSIM=0.2868

It can be seen from 1~table of table 2 at method filtering of the noisy medical image of different noise criteria differences using the present invention The PSNR and SSIM of the image obtained after reason are greatly increased both with respect to untreated noisy medical image, illustrate Model of the present invention has a good filtering performance, and filtering image processed by the invention can effective approaching to reality image.

Referring to Fig. 1~Fig. 8, it can be seen that noisy medical image is filtered using mean filter method and median filtering method The image visual effect that ripple processing is obtained is poor, illustrates that it still contains much noise, filtering performance is poor, filtered image has The obvious phenomenon that degrades, and use the method for the present invention to be filtered the image visual effect that processing is obtained to noisy medical image Have clear improvement, texture structure is effectively restored, while effectively inhibiting noise.The denoising method of the present invention is to different noises The medical image of standard deviation has excellent denoising effect.

Claims

1. a kind of Medical Image Denoising method based on resonance subspace analysis, it is characterised in that specific real according to following steps Apply：

Step 1：The medical image that a width contains random noise is inputted, this noisy medical image is designated as Y, Y ∈ R^M×N, Y=X+E, X For original medical image, E is the white Gaussian noise of addition, and M × N represents the size of image；

2.4, calculate the weighted value of the image block in similar block group；

2.5, in units of similar block group, initial filter, specific method are carried out to the target image block in each similar block group respectively It is as follows：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>{</mo> <msub> <mi>U</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>}</mo> <mo>=</mo> <mi>arg</mi> <mi> </mi> <mi>min</mi> <mfrac> <mn>1</mn> <mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <mo>|</mo> <mo>|</mo> <msub> <mi>P</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>U</mi> <mn>1</mn> </msub> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <msup> <msub> <mi>V</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&eta;</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>f</mi> <mn>1</mn> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>&gamma;</mi> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <msub> <mi>S</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>i</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>-</mo> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>j</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced>

Constraints

Wherein, η, γ are preset constant, K₁Represent the image number of blocks included in similar block group, P_kRepresent in similar block group k-th Image block,H is denoising parameter, U₁={ u₁,u₂,...u_d, U₁∈R^n×d, V₁={ v₁, v₂,...v_d, V₁∈R^n×d, U₁,V₁For left and right basic matrix, it is made up of orthogonal base vectors, f₁Estimate for the filtering of image block,Generation The weighted value of each image block in table similar block group,Square of F norms is represented, | | | |_TVRepresent TV norms；

2.5.4, the filter result of each target image block is added to the one and noisy medical image Y blank matrix phase with size On the position answered, the number average value obtained after each filtering for taking each pixel is obtained as the last estimate of this pixel To the initial filter result to noisy medical image；

3.1, according to the method for step 2.1~step 2.4 successively by last time filtered noisy medical image piecemeal, estimate that it is made an uproar Sound standard deviation, builds similar block group, calculates the weighted value of each image block in similar block group；

3.2, in units of similar block group, the target image block in each similar block group is filtered again respectively, specific method It is as follows：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>{</mo> <msub> <mi>U</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>}</mo> <mo>=</mo> <mi>arg</mi> <mi> </mi> <mi>min</mi> <mfrac> <mn>1</mn> <mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <mo>|</mo> <mo>|</mo> <msub> <mi>M</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>U</mi> <mn>2</mn> </msub> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <msup> <msub> <mi>V</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&eta;</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>f</mi> <mn>2</mn> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>&gamma;</mi> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <msubsup> <mi>S</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> <mo>&prime;</mo> </msubsup> <mo>|</mo> <mo>|</mo> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>i</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>-</mo> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>j</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced>

Constraints

Wherein, K₂Represent the image number of blocks included in similar block group, M_kRepresent k-th of image block, M in similar block group_k∈R^n×n, U₂={ u '₁,u'₂,...u'_d},U₂∈R^n×d, V₂={ v '₁,v'₂,...v'_d},V₂∈R^n×d, U₂,V₂For left and right basic matrix, by just Base vector is handed over to constitute,H is denoising parameter, f₂Estimate for the filtering of image block,Represent image The matrix of the composition reciprocal of pixel filter times in block,J₀Represent using last filtered image as The image block that basis is divided, that is, the pending target image block of this filtering,Represent each figure in similar block group As the weighted value of block；

3.2.2, with alternating iteration multiplier method, the denoising model in solution procedure 3.2.1 corresponds to the U of each image block₂、V₂；

3.2.4, the filter result of each target image block is added to the one and noisy medical image Y blank matrix phase with size On the position answered, the number average value obtained after each filtering for taking each pixel is obtained as the last estimate of this pixel To the filter result again to noisy medical image；

Repeat repeatedly to filter image using the filtering method of step 3, filtering each time is obtained after filtering with the last time Image based on, each image block or so basic matrix is calculated to being filtered to image, until left and right basic matrix is obtained to convergence Medical image figure after final denoising.

2. a kind of Medical Image Denoising method based on resonance subspace analysis according to claim 1, it is characterised in that The calculation formula of noise criteria difference is as follows in step 2：

σ₁=median (W_HH)/0.6745；

Wherein, median (W_HH) represent that image decomposes through 2-d wavelet after, the wavelet coefficient of obtained high pass-high pass subspace.

3. a kind of Medical Image Denoising method based on resonance subspace analysis according to claim 1, it is characterised in that The calculation formula of the weighted value of each image block in step 2 in similar block group is as follows：

<mrow> <msubsup> <mi>w</mi> <mi>k</mi> <mn>1</mn> </msubsup> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>P</mi> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> <mi>&alpha;</mi> </msubsup> </mrow> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <mo>|</mo> <mo>|</mo> <msub> <mi>P</mi> <mi>i</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> <mi>&alpha;</mi> </msubsup> </mrow> </mfrac> <mo>,</mo> </mrow>

Wherein, α represents time number formulary in above formula.

4. a kind of Medical Image Denoising method based on resonance subspace analysis according to claim 1, it is characterised in that The calculation formula of the weighted value of each image block in step 3 in similar block group is as follows：

<mrow> <msubsup> <mi>w</mi> <mi>k</mi> <mn>2</mn> </msubsup> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>M</mi> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> <mi>&alpha;</mi> </msubsup> </mrow> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <mo>|</mo> <msub> <mi>M</mi> <mi>i</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> <mi>&alpha;</mi> </msubsup> </mrow> </mfrac> <mo>,</mo> </mrow>

Wherein, α represents time number formulary in above formula.

5. a kind of Medical Image Denoising method based on resonance subspace analysis according to claim 1, it is characterised in that The method that similar block group is built in step 2 is as follows：

B, then sets distance threshold τ_d=3 σ₁ ²n², it is sized successively centered on each pixel in target image block as x × x Search window N (s), calculates the other image blocks and target image block P in the range of each search window₀Euclidean distance, calculation formula is such as Under：

<mrow> <mi>d</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>Y</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>Y</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>,</mo> <mi>t</mi> <mo>&Element;</mo> <mi>N</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

D, similar block group is built into by the similar image block selectedP_k∈R^n×n, k=0,1,2 ..., K₁。

6. a kind of Medical Image Denoising method based on resonance subspace analysis according to claim 1, it is characterised in that U is solved in step 2.5.2₁,V₁Method it is as follows：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>{</mo> <msub> <mi>U</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>}</mo> <mo>=</mo> <mi>arg</mi> <mi> </mi> <mi>min</mi> <mfrac> <mn>1</mn> <mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <mo>|</mo> <mo>|</mo> <msub> <mi>P</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>U</mi> <mn>1</mn> </msub> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <msup> <msub> <mi>V</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&eta;</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>f</mi> <mn>1</mn> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>&gamma;</mi> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <msub> <mi>S</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>i</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>-</mo> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>j</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mo><</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>f</mi> <mn>1</mn> </msub> <mo>-</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <msubsup> <mi>w</mi> <mi>k</mi> <mn>1</mn> </msubsup> <msub> <mi>U</mi> <mn>1</mn> </msub> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <msup> <msub> <mi>V</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <mo>></mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mfrac> <mi>&mu;</mi> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <msub> <mi>f</mi> <mn>1</mn> </msub> <mo>-</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <msubsup> <mi>w</mi> <mi>k</mi> <mn>1</mn> </msubsup> <msub> <mi>U</mi> <mn>1</mn> </msub> <msup> <msub> <mi>U</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>1</mn> </msub> <msup> <msub> <mi>V</mi> <mn>1</mn> </msub> <mi>T</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced>

2.5.2.2, each variable of initial filter model is initialized, U is initialized using singular value decomposition SVD₁、V₁ForV₁ ¹, λ₁'s It is initialized as the null matrix with image block formed objectsf₁Target image block is initialized as, f is designated as₁ ¹；

(1) f is solved₁ ^g+1：

Solve f₁Method it is as follows：

Gradient descent method updates f₁Formula be：Z represents that gradient descent method updates f₁Iterations, Dt represents to spread step-length, f₁₍₁₎=f₁ ^g；

(2) alternating iteration Optimization SolutionV₁ ^g+1：

1. U is initialized₁,V₁For orthogonal matrix, it is designated as

2. solveM is alternating iteration Optimization SolutionV₁ ^g+1Iterations：

Fixed V₁, solve U₁Method it is as follows：

To formulaEigenvalues Decomposition, by spy The descending sequence of value indicative, the corresponding characteristic vector of d characteristic value, that is, constitute U before taking₁；

Willf₁=f₁ ^g+1,Substitute into, obtain

3. solve

Fixed U₁, solve V₁Method it is as follows：

To formulaCarry out Eigenvalues Decomposition, By the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value, as V before taking₁；

Willf₁=f₁ ^g+1,Substitute into, obtain

Wherein,L₁ =D₁-S₁,I is unit battle array,I=1,2 ..., K₁：

<mrow> <msub> <mi>O</mi> <msub> <mi>d</mi> <mn>1</mn> </msub> </msub> <mo>=</mo> <mo>&lsqb;</mo> <mn>0</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mn>0</mn> <mo>;</mo> <mn>...</mn> <mo>;</mo> <mn>0</mn> <mo>,</mo> <mn>..</mn> <mo>,</mo> <mn>0</mn> <mo>&rsqb;</mo> <mo>,</mo> <msub> <mi>O</mi> <msub> <mi>d</mi> <mn>1</mn> </msub> </msub> <mo>&Element;</mo> <msup> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mn>1</mn> </msub> <mo>&times;</mo> <msub> <mi>d</mi> <mn>1</mn> </msub> </mrow> </msup> <mo>;</mo> </mrow>

4., with For next iteration process U₁,V₁Initial value, repeat step 2., 3., until convergence, obtainV₁ ^g+1；

<mrow> <msubsup> <mi>&lambda;</mi> <mn>1</mn> <mrow> <mi>g</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>&lambda;</mi> <mn>1</mn> <mi>g</mi> </msubsup> <mo>+</mo> <mi>&mu;</mi> <mrow> <mo>(</mo> <msubsup> <mi>f</mi> <mn>1</mn> <mrow> <mi>g</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>1</mn> </msub> </munderover> <msubsup> <mi>w</mi> <mi>k</mi> <mn>1</mn> </msubsup> <msubsup> <mi>U</mi> <mn>1</mn> <mrow> <mi>g</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <msubsup> <mi>U</mi> <mn>1</mn> <mrow> <mi>g</mi> <mo>+</mo> <msup> <mn>1</mn> <mi>T</mi> </msup> </mrow> </msubsup> <msub> <mi>P</mi> <mi>k</mi> </msub> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mi>g</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mi>g</mi> <mo>+</mo> <msup> <mn>1</mn> <mi>T</mi> </msup> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow>

(4) withV₁ ^g+1、f₁ ^g+1、For next iteration process U₁、V₁、f₁、λ₁Initial value, repeat step (1)~(3), Until the convergence of each variable, obtains final U₁,V₁。

7. a kind of Medical Image Denoising method based on resonance subspace analysis according to claim 1, it is characterised in that U is solved in step 3.2.2₂、V₂Method it is as follows：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>{</mo> <msub> <mi>U</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>}</mo> <mo>=</mo> <mi>arg</mi> <mi> </mi> <mi>min</mi> <mfrac> <mn>1</mn> <mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <mo>|</mo> <mo>|</mo> <msub> <mi>M</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>U</mi> <mn>2</mn> </msub> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <msup> <msub> <mi>V</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&eta;</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>f</mi> <mn>2</mn> </msub> <mo>|</mo> <msub> <mo>|</mo> <mrow> <mi>T</mi> <mi>V</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>&gamma;</mi> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <msubsup> <mi>S</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> <mo>&prime;</mo> </msubsup> <mo>|</mo> <mo>|</mo> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>i</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>-</mo> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>j</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mo><</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>f</mi> <mn>2</mn> </msub> <mo>-</mo> <mi>A</mi> <mo>-</mo> <mover> <mi>b</mi> <mo>~</mo> </mover> <mo>.</mo> <mo>*</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <msubsup> <mi>w</mi> <mi>k</mi> <mn>2</mn> </msubsup> <msub> <mi>U</mi> <mn>2</mn> </msub> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <msup> <msub> <mi>V</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <mo>></mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mfrac> <mi>&mu;</mi> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <msub> <mi>f</mi> <mn>2</mn> </msub> <mo>-</mo> <mi>A</mi> <mo>-</mo> <mover> <mi>b</mi> <mo>~</mo> </mover> <mo>.</mo> <mo>*</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <msubsup> <mi>w</mi> <mi>k</mi> <mn>2</mn> </msubsup> <msub> <mi>U</mi> <mn>2</mn> </msub> <msup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <msub> <mi>M</mi> <mi>k</mi> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <msup> <msub> <mi>V</mi> <mn>2</mn> </msub> <mi>T</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mi>F</mi> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced>

Wherein, μ is default constant；

3.2.2.2, each variable of denoising model in initialization step 3.2, U is initialized using singular value decomposition SVD₂,V₂ForBy λ₂It is initialized as the null matrix with image block formed objectsf₂It is initialized as the phase of last filtered image Like the target image block in block group, it is designated as

3.2.2.3.1, solve

Solve f₂Method it is as follows：

Gradient descent method updates f₂Formula be：L represents that gradient descent method updates f₂Iterations, Dt represents diffusion step-length,

Simultaneous above formula：WillSubstitution formula is produced

3.2.2.3.2, solve

(1) U is initialized₂,V₂For orthogonal matrix, it is designated asξ be initialized as with tile size identical null matrix, be designated as ξ¹, alternately solve t, V₂、U₂, ξ iterations be designated as d；

(2) t is sought^d：

U is fixed using least square method₂、V₂Solve t method be：

By expression formulaCalculated using least square method and obtain T, then this vector is changed For the matrix of corresponding form, t is produced；

Wherein, F, T, R, b₂Respectively matrixVectorization form, B_bFor diagonal matrix, element is vector b on its diagonal₂The element of middle correspondence position, b₁Corresponded to for target image block in B matrixes The matrix of position, b₁∈R^n×n, B matrixes are the matrix of the corresponding filter times composition of entire image each pixel；

Willξ=ξ^dSubstitute into, produce t^d；

(3) solve

Solve first

B, is solvedP is alternative method solution in step 3.2.2.3.2Iterations：

1. solve

Fixed V₂Seek U₂Method it is as follows：

To formulaCarry out characteristic value point Solution, by the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value before taking constitutes U₂；

By t=t^d,ξ=ξ^dSubstitute into, obtain

2. solve

Fixed U₂Seek V₂Method it is as follows：

To formulaCarry out characteristic value Decompose, by the descending sequence of characteristic value, the corresponding characteristic vector of d characteristic value before taking constitutes V₂；

By t=t^d,ξ=ξ^dSubstitute into, obtain

Wherein, ξ is Lagrange multiplier, and κ is preset constant, <mrow> <msub> <mi>Y</mi> <msub> <mi>V</mi> <mn>2</mn> </msub> </msub> <mo>=</mo> <mo>&lsqb;</mo> <msub> <mi>M</mi> <mn>0</mn> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>M</mi> <mn>1</mn> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msub> <mi>M</mi> <msub> <mi>K</mi> <mn>2</mn> </msub> </msub> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>&rsqb;</mo> <mo>,</mo> <msub> <mi>Y</mi> <msub> <mi>U</mi> <mn>2</mn> </msub> </msub> <mo>=</mo> <mo>&lsqb;</mo> <msubsup> <mi>M</mi> <mn>0</mn> <mi>T</mi> </msubsup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mo>,</mo> <msubsup> <mi>M</mi> <mn>1</mn> <mi>T</mi> </msubsup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msubsup> <mi>M</mi> <msub> <mi>K</mi> <mn>2</mn> </msub> <mi>T</mi> </msubsup> <msub> <mi>U</mi> <mn>2</mn> </msub> <mo>&rsqb;</mo> <mo>,</mo> </mrow> L2 =D2-S2, I are unit battle array, i=1,2 ..., K2：

<mrow> <msub> <mi>O</mi> <msub> <mi>d</mi> <mn>2</mn> </msub> </msub> <mo>=</mo> <mo>&lsqb;</mo> <mn>0</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mn>0</mn> <mo>;</mo> <mn>...</mn> <mo>;</mo> <mn>0</mn> <mo>,</mo> <mn>..</mn> <mo>,</mo> <mn>0</mn> <mo>&rsqb;</mo> <mo>,</mo> <msub> <mi>O</mi> <msub> <mi>d</mi> <mn>2</mn> </msub> </msub> <mo>&Element;</mo> <msup> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mn>2</mn> </msub> <mo>&times;</mo> <msub> <mi>d</mi> <mn>2</mn> </msub> </mrow> </msup> <mo>;</mo> </mrow>

C, withFor the U of next iteration process₂、V₂Initial value, repeat step 1., 2. until convergence, Obtain

(4) ξ is solved^d+1：

<mrow> <msup> <mi>&xi;</mi> <mrow> <mi>d</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <msup> <mi>&xi;</mi> <mi>d</mi> </msup> <mo>+</mo> <mi>&kappa;</mi> <mrow> <mo>(</mo> <msub> <mi>b</mi> <mn>1</mn> </msub> <mo>.</mo> <mo>*</mo> <msup> <mi>t</mi> <mi>d</mi> </msup> <mo>-</mo> <msubsup> <mi>U</mi> <mrow> <mn>2</mn> <mrow> <mo>(</mo> <mi>d</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <msup> <msubsup> <mi>U</mi> <mrow> <mn>2</mn> <mrow> <mo>(</mo> <mi>d</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mi>T</mi> </msup> <msubsup> <mi>EV</mi> <mrow> <mn>2</mn> <mrow> <mo>(</mo> <mi>d</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <msup> <msubsup> <mi>V</mi> <mrow> <mn>2</mn> <mrow> <mo>(</mo> <mi>d</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mi>T</mi> </msup> <mo>)</mo> </mrow> </mrow>

<mrow> <msubsup> <mi>&lambda;</mi> <mn>2</mn> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>&lambda;</mi> <mn>2</mn> <mi>c</mi> </msubsup> <mo>+</mo> <mi>&mu;</mi> <mrow> <mo>(</mo> <msubsup> <mi>f</mi> <mn>2</mn> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <mi>A</mi> <mo>-</mo> <mover> <mi>b</mi> <mo>~</mo> </mover> <mo>.</mo> <mo>*</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <msub> <mi>K</mi> <mn>2</mn> </msub> </munderover> <msubsup> <mi>w</mi> <mi>k</mi> <mn>2</mn> </msubsup> <msubsup> <mi>U</mi> <mn>2</mn> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <msubsup> <mi>U</mi> <mn>2</mn> <mrow> <mi>c</mi> <mo>+</mo> <msup> <mn>1</mn> <mi>T</mi> </msup> </mrow> </msubsup> <msub> <mi>M</mi> <mi>k</mi> </msub> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mi>c</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mi>c</mi> <mo>+</mo> <msup> <mn>1</mn> <mi>T</mi> </msup> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow> 6

3.2.2.3.4, calculate what is obtained with above-mentioned stepsFor the U during next iteration₂、 V₂、f₂、λ₂Initial value, repeat step 3.2.2.3.1~3.2.2.3.3 until each variable restrain, obtain final U₂,V₂。