CN100387194C - Interation curative wave filtration combined weighted least squares positron emission tomography method - Google Patents

Abstract

The present invention discloses an iteration curvature flow filtering weighted least square positron emission tomography method. The present invention comprises the following steps: firstly getting observation data; determining the scale of the original image; comparing each component in the observation data with 1; dividing each component in the error variance and filtering an intermediate image with curvature flow and getting a filtered image as the original image; repeating the process until a reconstructed image is contracted and getting the finished image. The curvature flow filtering is added in the least square iterating method by the characteristics of the features of the current PET image-forming data, which inhibits pattern noise and improves the image quality; the present invention has the advantages of simple realization, high treating speed and good effect on inhibiting the background interference of the reconstructed image.

Description

The interation curative wave filtration combined weighted least squares positron emission tomography method

Technical field

The invention belongs to positron emission tomography (PET:Positron emission tomography) technical field of imaging, the invention particularly relates to a kind of interation curative wave filtration combined weighted least squares positron emission tomography method that is used to make up image.

Background technology

In modern PET imaging system, the precorrection that meets at random of original Poisson observed data is operated, especially postpone the window coincidence correction, become an essential important step.It can suppress the major part signal to noise ratio of coincident event (especially ambient interferences) raising observed data at random effectively.Yet because the processing of precorrection makes observed data not obey original Poisson statistics characteristic.Like this based on the method for estimation of original Poisson maximum likelihood not feasible.In this case, we consider weighted least require method usually.

The observed data model of existing least square method is:

$y_i=\underset{j}{\mathrm{\Σ}}{a}_{\mathrm{ij}}{f}_i+{\mathrm{\ϵ}}_i---\left(1\right)$

Y wherein _iI the number of photons that probe access detects of expression precorrection, 1≤i≤M, M are total probe access number; f _jRepresent the number of photons that j pixel goes out to send, f _j〉=0,1≤j≤N, N are total number of picture elements; a _IjRepresent that photon energy that j pixel go out to send is by i the detected probability of probe access; ε _iExpression is attached to the random error on the photon number that i probe access detect, its statistical property obedience zero-mean, and variance is d _iNormal distribution.According to (1), the object function that we rebuild can be expressed as:

$J\left(f\right)=\arg \underset{f}{\min }{(y-\mathrm{Af})}^T{\mathrm{\Σ}}^{-1}(y-\mathrm{Af}),$ And Condition of Non-Negative Constrains: f 〉=0 (2)

Here y is by y _i, the observed data vector that 1≤i≤M forms; F is by f _j, the image vector that 1≤j≤N forms; A is by a _IjSystem's probability matrix of the M * N that forms; ∑ is M * M weighting diagonal matrix, and i element on its diagonal is d _i

Choosing of weighting matrix is extremely important, and its quality directly has influence on quality of reconstructed images.We can adopt following variance method of estimation:

d _i＝max(y _i，1)，1≤i≤M，(3)

Under Condition of Non-Negative Constrains, if image

Separate for the weighted least-squares of object function (2), must satisfy the Kuhn-Tucker condition so:

a) ${\hat{f}}_j\underset{i}{\mathrm{\Σ}}\frac{1}{d_i}{a}_{\mathrm{ij}}(y_i-\underset{p}{\mathrm{\Σ}}{a}_{\mathrm{ip}}{\hat{f}}_p)={\hat{f}}_j\underset{i}{\mathrm{\Σ}}\frac{1}{d_i}{a}_{\mathrm{ij}}{y}_i,$ If ${\hat{f}}_j>0,j=1,...,N---\left(4\right)$

b) ${\hat{f}}_j\underset{i}{\mathrm{\Σ}}\frac{1}{d_i}{a}_{\mathrm{ij}}(y_i-\underset{p}{\mathrm{\Σ}}{a}_{\mathrm{ip}}{\hat{f}}_p)\≤{\hat{f}}_j\underset{i}{\mathrm{\Σ}}\frac{1}{d_i}{a}_{\mathrm{ij}}{y}_i,$ If ${\hat{f}}_j\≤0,j=1,...,N---\left(5\right)$

We can directly obtain existing least square alternative manner according to (4):

${\hat{f}}_j^{n+1}={\hat{f}}_j^n\frac{{\mathrm{\Σ}}_i{a}_{\mathrm{ij}}{y}_i/d_i}{{\mathrm{\Σ}}_i{a}_{\mathrm{ij}}{\mathrm{\Σ}}_p{a}_{\mathrm{ip}}{\hat{f}}_p^n/d_i},j=1,...,N,n=\mathrm{0,1,2},...\left(6\right)$

N represents iterations herein.

Because the incompleteness of data for projection itself and the ill-posedness of method for reconstructing only use weighted least require method to can not get ideal reconstructed image to a great extent.Phenomenons such as it is rough that its result is easy to generate image, and noise is more obvious.Common practices is the regularization technology, considers simultaneously that promptly the regular terms relevant with image smoothing suppresses effect of noise.Yet the complexity of its realization is different with the difference of regular terms.For example, some regularization methods can bring good reconstruction effect, but realize comparatively complexity, and level and smooth dynamics is difficult to control.The introducing of regular terms can not suppress background noise preferably in addition.For the shortcoming of prior art is described, we are according to the template data shown in Figure 1 emulation that experimentizes.Fig. 2 has provided and has utilized emulated data, according to the result of existing weighted least require method reconstruction.We can find out significantly that reconstructed image is not ideal enough, and making an uproar phenomenon is comparatively obvious.Therefore how reconstructing high-quality PET image is one of problem that presses at present solution.

Summary of the invention

The invention provides a kind of interation curative wave filtration combined weighted least squares positron emission tomography method that can improve reconstructed image quality.

The present invention adopts following technical scheme:

A kind of interation curative wave filtration combined weighted least squares positron emission tomography method that is used to make up image, adopt the following step:

1) obtains observed data, write down its scale M, and it is kept at vectorial y=[y ₁..., y _M] ^TIn; According to the dimensional requirement for the treatment of reconstructed image, determine the scale N of initial pictures, given its initial gray value is greater than 1, and it is kept at vector f=[f ₁..., f _N] ^TIn,

2) according to the size of observed data y and initial pictures f, calculate the probability matrix A of system of M * N,

3) each component among the observed data y is carried out one by one with 1 relatively: if than 1 greatly then keep, then be set to 1 for a short time than 1, and comparative result be kept among the vectorial d d=[d ₁..., d _M] ^T, as the error variance of observed data,

4) with each component among the observed data y and the corresponding one by one observation data for projection that is divided by and obtains revising of each component among the error variance d

$\tilde{y}={[{\tilde{y}}_1,...,{\tilde{y}}_M]}^T,$ Wherein ${\tilde{y}}_i=y_i/d_i,$ i＝1，…，M，

5) with the transposition of the probability matrix A of system and the observed data of correction Multiplying each other obtains vectorial g, $g=A^T\tilde{y},$ g＝[g ₁，…，g _N] ^T，

6) probability matrix A of system and initial pictures f are multiplied each other, obtain forward projection p, p=[p ₁... p _M] ^T, again with each component among this forward projection p and the corresponding one by one forward projection that is divided by and obtains revising of each component among the error variance d

$\tilde{p}={[{\tilde{p}}_1,...,{\tilde{p}}_M]}^T,$ Wherein ${\tilde{p}}_i=p_i/d_i,$ I=1 ..., M utilizes transposition and this forward projection of the probability matrix A of system then

Multiply each other and obtain the h of back projection: $h=A^T\tilde{p},$ h＝[h ₁，…，h _N] ^T，

7) each component of vectorial g and each component corresponding one by one being divided by among the h of back projection, obtained image correction value c:c=[c ₁..., c _N] ^T, c wherein _j=g _j/ h _j, j=1 ..., N,

8) each component among this image correction value c is obtained intermediate image with corresponding one by one the multiplying each other of each component among the initial pictures f

$\hat{f}={[{\hat{f}}_1,...,{\hat{f}}_N]}^T,$ Wherein ${\hat{f}}_j=c_j\×f_j,$ j=1，…，N，

9) with intermediate image

Carry out curvature flow filter at least one time, can obtain filtering image

10) with this filtering image

As initial pictures, turned back to for the 6th step, repeat the image convergence of this process after rebuilding, can obtain the image of final imaging.

The present invention is in conjunction with the characteristics of current PE T imaging data, and the drawback at existing method exists has proposed a kind of new solution: the interation curative wave filtration combined weighted least squares formation method.The present invention adds curvature flow filter mechanism in existing least square method iterative process, suppress making an uproar of image with this, improves the quality of image.

Compared with prior art, the outstanding effect that the present invention has is: realize comparatively simply, wherein the curvature flow filter process acts directly on the two dimensional image, needn't carry out conversion, so processing speed is very fast; The number of times that can regulate filtering to the control of the smoothness of required reconstructed image is accomplished.Because the filtering of curvature flow is a kind of preferable image filtering technique, it does not influence the original basic feature of image to smoothed image the time, such as the edge, and line, angle etc., therefore the image of rebuilding can keep due basic feature information in the original image equally; In addition, the ambient interferences that suppresses reconstructed image also there is good effect.

Fig. 3 has provided the result that the present invention rebuilds emulated data, wherein carries out Filtering Processing between 3 iteration.Itself and Fig. 2 are compared, and we can clearly find, the clear picture of rebuilding by method of the present invention, smooth can reflect the original feature of original image.

For effectiveness of the present invention there being a quantitative recognition, we weigh the quality of this formation method with the reconstructed image signal-to-noise ratio (SNR) Criterion.This criterion is defined as follows:

$\mathrm{SNR}\left(\mathrm{dB}\right)={10\log }_{10}\left(\frac{{||f||}^2}{{||f^{*}-f||}^2}\right)---\left(7\right)$

F wherein ^*Expression simulation modular view data.It is just high more to rebuild signal to noise ratio, illustrates that then formation method is good more.By calculating, we list the reconstruction signal to noise ratio of Fig. 2 and Fig. 3 respectively in the table 1, to make comparisons.

Table 1

Obviously, can obtain higher signal to noise ratio with the method for the invention, its value exceeds 4dB nearly than existing weighted least require method.This has illustrated significantly that it is superior to the general property taken advantage of weighted least require method.

Table 2

The present invention rebuilds signal to noise ratio to be increased and increases along with the number of times of filtering.But the amplitude that increases reduces gradually along with the increase of number of times, so filter times can control in certain scope, such as 10 times.Show that in addition the method that we invent converges to consistent reconstructed results the most at last.

In order further to verify the superiority of the method that we invent, we are illustrated with actual 3 clinical datas that obtain.Fig. 4, Fig. 5, Fig. 6 use the prior art reconstructed results.Clearly, picture quality is relatively poor, is not easy to clinical research.Fig. 7, Fig. 8 and Fig. 9 be corresponding diagram 4, Fig. 5 and Fig. 6 respectively, and they are the results with method imaging of the present invention.Therefrom we can clearly find, interation curative wave filtration combined weighted least squares method image quality of the present invention is quite superior, and it is not only level and smooth image has kept the edge, and also quite effective to suppressing ambient interferences.Therefore the present invention is more suitable for the clinical diagnosis needs.

Description of drawings

Fig. 1 is the thoracic cavity template image that is used for testing formation method

The result of Fig. 2 for rebuilding with existing weighted least require method

Fig. 3 is with the result's (3 filtering) after the imaging of the present invention

The result of Fig. 4 for clinical data 1 being carried out imaging with existing weighted least require method

The result of Fig. 5 for clinical data 2 being carried out imaging with existing weighted least require method

The result of Fig. 6 for clinical data 3 being carried out imaging with existing weighted least require method

The result (5 filtering) of Fig. 7 for clinical data 1 being carried out imaging with the present invention

The result (5 filtering) of Fig. 8 for clinical data 2 being carried out imaging with the present invention

The result (5 filtering) of Fig. 9 for clinical data 3 being carried out imaging with the present invention

The specific embodiment

The specific embodiment of the present invention is as follows:

1) obtains the observed data of carrying out coincidence correction at random, determine its scale M, and save as at vectorial y=[y ₁..., y _M] ^TAccording to the dimensional requirement for the treatment of reconstructed image, determine the scale N of initial pictures, given its initial gray value is greater than 1, and it is kept at vector f=[f ₁..., f _N] ^TIn.

2), calculate the probability matrix A of system of M * N according to the size of observed data y and initial pictures f.In order to calculate this system's probability matrix, we answer the geometry of known scanner, then can be according to certain method, and for example visual angle method (Angle of view) is calculated.

In practical situation, the observed data collection is from PET video picture scanner.For the quality of appraisal procedure, we also can generate the observed data that needs by computer simulation experiment.Here we illustrate with the PET thoracic cavity template of a Computer Simulation shown in Figure 1.This thoracic cavity template size is 128 * 128 PEL matrix, i.e. N=128 * 128.We suppose that the data for projection scale is M=192 * 192, and it represents 192 projection angles, and 192 probe access are arranged on each angle.We suppose that further passage and interchannel spacing are a pixel, are convenient to our computing system probability matrix A like this.

The method of coincidence correction is a lot of at random, and modern PET scanner all adopts and postpones the window coincidence correction.By its principle, we can make following emulation experiment to generate required observed data.We multiply each other with system's probability matrix and thoracic cavity template image earlier and obtain muting observed data, then all data summations are obtained total number of photons.Total number of photons is counted M divided by total probe access, estimate the average photon number that each probe access is caught.Get 25% random mark average (the 25%th, according to just before giving birth the empirical value of the analysis gained of data) as a setting of this meansigma methods.We add each noiseless observed data the average of the Poisson observation data for projection that this background random mark average is used as emulation.In conjunction with the Poisson PRNG on the computer, produce required uncorrected observation data for projection with this average.With same background random mark average, in conjunction with PRNG, we can simulate other one group of Poisson data and meet data at random as postponing window.Uncorrected observation data for projection is deducted the observation data for projection y that these delay window random mark data obtain final precorrection.

3) utilize observation data for projection y, the error variance d=[d of the method calculating observation data of being given in conjunction with formula (3) ₁..., d _M] ^T

4) with each component among the observed data y and the corresponding one by one observation data for projection that is divided by and obtains revising of each component among the error variance d

$\tilde{y}=[{{\tilde{y}}_1,...,{\tilde{y}}_M]}^T;$ Wherein ${\tilde{y}}_i=y_i/d_i,$ i＝1，…，M，

5) with the transposition of the probability matrix A of system and the observed data of correction

Multiplying each other obtains vectorial g, that is: $g=A^T\tilde{y},$ g＝[g ₁，…，g _N] ^T，

6) probability matrix A of system and initial pictures f are multiplied each other, obtain forward projection p, p=[p ₁..., p _M] ^T, again with each component among this forward projection p and the corresponding one by one forward projection that is divided by and obtains revising of each component among the error variance d

$\tilde{p}={[{\tilde{p}}_1,...{,\tilde{p}}_M]}^T,$ Wherein ${\tilde{p}}_i=p_i/d_i,$ i＝1，…，M。Utilize transposition and this forward projection of the probability matrix A of system then

Multiply each other and obtain the h of back projection: $h=A^T\tilde{p},$ h＝[h ₁，…，h _N] ^T，

7) each component of vectorial g and each component corresponding one by one being divided by among the h of back projection, obtained image correction value c:c=[c ₁..., c _N] ^T, c wherein _j=g _j/ h _j, j=1 ..., N,

8) each component among this image correction value c is obtained intermediate image with corresponding one by one the multiplying each other of each component among the initial pictures f

$\hat{f}={[{\hat{f}}_1,...,{\hat{f}}_N]}^T,$ Wherein ${\hat{f}}_j=c_j\×f_j,$ j＝1，…，N，

9) middle image is carried out curvature flow filter: curvature flow filter can effectively be removed noise, but does not damage image original characteristics information, for example edge, line, angle etc.Its filtering can be represented with following partial differential equation:

$\frac{\∂f(h,v;t)}{\∂t}=\▿\·\left(\frac{\▿f(h,v;t)}{|\▿f(h,v;t)|}\right)---\left(7\right)$

And boundary condition: $\frac{\∂f(h,v;t)}{\∂\stackrel{\→}{N}}{|}_{f=\∂\mathrm{\Ω}}=0,$ Initial condition: f (h, v; 0)=f ⁰Here

The expression gradient operator; For the ease of difference, we use h, the level of v presentation video and the coordinate of vertical both direction; Ω is the image domain of definition;

Be image boundary,

Be the boundary method vector; The t express time.Introducing time evolution step delta t (it is comparatively suitable to get Δ t=0.0005 here) also launches equation (7) the right item:

$f^{t+1}\text{=}{f}^t+\mathrm{\Δt}\×\left(\frac{f_{\mathrm{hh}}^t{\left(f_v^t\right)}^2-2f_{\mathrm{hv}}^t{f}_h^t{f}_v^t+f_{\mathrm{vv}}^t{\left(f_h^t\right)}^2}{\sqrt{{\left(f_h^t\right)}^2+{\left(f_v^t\right)}^2}}\right)---\left(8\right)$

F wherein ^t≡ f (h, v; T); $f_h^t=\∂f(h,v;t)/\∂h,$ $f_v^t\≡\∂f(h,v;t)/\∂v$ Be the level of t moment function f and the single order partial derivative of vertical both direction; $f_{\mathrm{hh}}^t\≡{\∂}^{2}f(h,v;t)/{\∂h}^2,$ $f_{\mathrm{hv}}^t\≡\∂f(h,v;t)/\∂h\∂v$ With $f_{\mathrm{vv}}^t\≡{\∂}^{2}f(h,v;t)/{\∂v}^2$ Be second-order partial differential coefficient.

In practical situation, we consider discrete picture usually.For using formula (8) carries out filtering, we need introduce suitable difference scheme and come partial derivative in the replacement formula (8).Narration for convenience, we can suppose that image is that (square formation and the pixel pitch of K * K=N) are unit length 1 to K * K.Note (f _h ^t) _jBe function f (h, v; T) at point (h _j, v _j) (1≤h wherein _j≤ K, 1≤v _jThe single order partial derivative of≤h direction K) located, then it can be represented with following first-order difference:

${\left(f_h^t\right)}_j=f(h_j+1,v_j;t)-f(h_j,v_j;t),j=1,...,N---\left(9\right)$

Wherein: h _j=int[j/K]+1, v _j=j-int[j/N] * N (int[x] the expression integer part of getting x) and j=(h _j-1) * K+v _jIn the same way, we can also obtain:

${\left(f_v^t\right)}_j=f(h_j,v_j+1;t)-f(h_j,v_j;t),---\left(10\right)$

${\left(f_{\mathrm{hh}}^t\right)}_j=f(h_j+1,v_j;t)-2f(h_j,v_j;t)+f(h_j-1,v_j;t),---\left(11\right)$

${\left(f_{\mathrm{vv}}^t\right)}_j=f(h_j,v_j+1;t)-2f(h_j,v_j;t)+f(h_j,v_j-1;t),---\left(12\right)$

${\left(f_{\mathrm{hv}}^t\right)}_j=\frac{(f(h_j+1,v_j+1;t)+f(h_j-1,v_j-1;t))-(f(h_j-1,v_j+1;t)+f(h_j+1,v_j-1;t))}{4}---\left(13\right)$

Order in like manner $f_j^t=f(h_j,v_j;t),$ j＝1，…，N。We can be rewritten as (8) like this:

${f_j}^{t+1}={f_j}^i+\mathrm{\Δt}\×\frac{{\left(f_{\mathrm{hh}}^t\right)}_j{\left(f_v^t\right)}_j^2-2{\left(f_{\mathrm{hv}}^t\right)}_j{\left(f_h^t\right)}_j{\left(f_v^t\right)}_j+{\left(f_{\mathrm{vv}}^t\right)}_j{\left(f_h^t\right)}_j^2}{\sqrt{{\left(f_h^t\right)}_j^2+{\left(f_v^t\right)}_j^2}},j=1,...,N,---\left(14\right)$

And t=0 ..., t _Max, t here _MaxFor total time is filter times.For the image border, we do the edge replication processes, promptly for formula (9) to formula (13):

If h _j-1≤0, then put h _j-1 is 1; If v _j-1≤0, then put v _j-1 is 1;

If h _j+ 1 〉=K then puts h _j+ 1 is K; If v _j+ 1 〉=K then puts v _j+ 1 is K.

According to above-mentioned curvature flow filter principle, to middle image

Carrying out curvature flow filter can followingly operate: total filter times t is set _Max(general 1≤t _Max≤ 10 is comparatively reasonable); Initial condition is set: $f_j^0={\hat{f}}_j,$ J=1 ..., N, and put t=0; By formula (14) carry out filtering, put t=t+1; Recurring formula (14) is until t=t _MaxFinish.Final filtering result is kept at vector

In, $\tilde{f}={[{\tilde{f}}_1,...,{\tilde{f}}_N]}^T.$

10) with this filtering image

As initial pictures f, turned back to for the 6th step, repeat the image convergence of this process after rebuilding, can obtain the image of final imaging.

Claims (1)

1. interation curative wave filtration combined weighted least squares positron emission tomography method that is used to make up image is characterized in that adopting the following step:

1) obtains observed data, write down its scale M, and it is kept at vectorial y=[y ₁..., y _M] ^TIn; According to the dimensional requirement for the treatment of reconstructed image, at square-shaped image, the length of side of establishing it is K, determines the scale N=K * K of initial pictures, and given its initial gray value is greater than 1, and it is kept at vector f=[f ₁..., f _N] ^TIn, the total degree I of setting iteration, I=3,

2) according to the size of observed data y and initial pictures f, calculate the probability matrix A of system of M * N,

3) each component among the observed data y is carried out one by one with 1 relatively: if than 1 greatly then keep, then be set to 1 for a short time than 1, and comparative result be kept among the vectorial d d=[d ₁..., d _M] ^T, as the error variance of observed data,

4) with each component among the observed data y and the corresponding one by one observation data for projection that is divided by and obtains revising of each component among the error variance d $\tilde{y}:\tilde{y}={[{\tilde{y}}_1,...,{\tilde{y}}_M]}^T,$ Wherein ${\tilde{y}}_i=y_i/d_i,$ i＝1，…，M，

5) with the transposition of the probability matrix A of system and the observed data of correction

Multiplying each other obtains vectorial g, $g=A^T\tilde{y},$ g＝[g ₁，…，g _N] ^T，

6) m=0 is set,

7) probability matrix A of system and initial pictures f are multiplied each other, obtain forward projection p, p=[p ₁..., p _M] ^T, again with each component among this forward projection p and the corresponding one by one forward projection that is divided by and obtains revising of each component among the error variance d

: $\tilde{p}={[{\tilde{p}}_1,...,{\tilde{p}}_M]}^T,$ Wherein ${\tilde{p}}_i=p_i/d_i,$ I=1 ..., M utilizes transposition and this forward projection of the probability matrix A of system then

Multiply each other and obtain the h of back projection: $h=A^T\tilde{p},$ h＝[h ₁，…，h _N] ^T，

8) each component of vectorial g and each component corresponding one by one being divided by among the h of back projection, obtained image correction value c:c=[c ₁..., c _N] ^T, c wherein _j=g _j/ h _j, j=1 ..., N,

9) each component among this image correction value c is obtained intermediate image with corresponding one by one the multiplying each other of each component among the initial pictures f $\hat{f}:\hat{f}={[{\hat{f}}_1,...,{\hat{f}}_N]}^T,$ Wherein ${\hat{f}}_j=c_j\×f_j,$ j＝1，…，N，

10) t=0 is set, and filtering total degree t _Max, 1≤t _Max≤ 10,

11) with intermediate image

Again note is made a new intermediate images vector ${\hat{f}}^t:{\hat{f}}^t={[{\hat{f}}_l^t,...,{\hat{f}}_N^t]}^T,$

12) be calculated as follows the derivative vector with this intermediate images:

${\hat{f}}_h^t={[{\left(f_h^t\right)}_1,...,{\left(f_h^t\right)}_N]}^T,$ ${\hat{f}}_v^t={[{\left(f_v^t\right)}_1,...,{\left(f_v^t\right)}_N]}^T,$

${\hat{f}}_{\mathrm{hh}}^t={[{\left(f_{\mathrm{hh}}^t\right)}_1,...,{\left(f_{\mathrm{hh}}^t\right)}_N]}^T,$ ${\hat{f}}_{\mathrm{vv}}^t={[{\left(f_{\mathrm{vv}}^t\right)}_1,...,{\left(f_{\mathrm{vv}}^t\right)}_N]}^T,$

${\hat{f}}_{\mathrm{hv}}^t={[{\left(f_{\mathrm{hv}}^t\right)}_1,...,{\left(f_{\mathrm{hv}}^t\right)}_N]}^T,$

Wherein

${\left(f_h^t\right)}_j={\hat{f}}_{h_j\×K+v_j}^t-{\hat{f}}_{(h_j-1)\×K+v_j}^t,$

${\left(f_v^t\right)}_j={\hat{f}}_{(h_j-1)\×K+v_j+1}^t-{\hat{f}}_{(h_j-1)\×K+v_j}^t,$

${\left(f_{\mathrm{hh}}^t\right)}_j={\hat{f}}_{h_j\×K+v_j}^t-2\×{\hat{f}}_{(h_j-1)\×K+v_j}^t+{\hat{f}}_{(h_j-2)\×K+v_j}^t,$

${\left(f_{\mathrm{vv}}^t\right)}_j={\hat{f}}_{(h_j-1)\×K+v_j+1}^t-2\×{\hat{f}}_{(h_j-1)\×K+v_j}^t+{\hat{f}}_{(h_j-1)\×K+v_j-1}^t,$

${\left(f_{\mathrm{hv}}^t\right)}_j=\frac{1}{4}\×({\hat{f}}_{h_j\×K+v_j+1}^t+{\hat{f}}_{(h_j-2)\×K+v_j-1}^t-{\hat{f}}_{(h_j-2)\×K+v_j+1}^t-{\hat{f}}_{h_j\×K+v_j-1}^t),$

J=1 ..., N, the h here _jAnd v _jCalculate with following two formulas:

h _j＝int[j/K]+1，v _j＝j-int[j/k]×K

Int[j/K herein] expression gets the merchant of j divided by K,

F wherein ^t≡ f (h, v; T); $f_h^t\≡\∂f(h,v;t)/\∂h,$ $f_v^t\≡\∂f(h,v;t)/\∂v$ Be the level of t moment function f and the single order partial derivative of vertical both direction; $f_{\mathrm{hh}}^t\≡{\∂}^{2}f(h,v;t)/{\∂h}^2,$ $f_{\mathrm{hv}}^t\≡\∂f(h,v;t)/\∂h\∂v$ With $f_{\mathrm{vv}}^t\≡{\∂}^{2}f(h,v;t)/{\∂v}^2$ Be second-order partial differential coefficient,

13) calculate middle filtering image ${\tilde{f}}^t:{\tilde{f}}^t={[{\tilde{f}}_1^t,...,{\tilde{f}}_N^t]}^T,$ Wherein

${\tilde{f}}_j^t={\hat{f}}_j^t+\mathrm{\Δt}\×\frac{{\left(f_{\mathrm{hh}}^t\right)}_j\×{\left((f_v^t{)}_j\right)}^2-2\×{\left(f_{\mathrm{hv}}^t\right)}_j\×{\left(f_h^t\right)}_j\×{\left(f_v^t\right)}_j+{\left(f_{\mathrm{vv}}^t\right)}_j\×{\left({\left(f_h^t\right)}_j\right)}^2}{\sqrt{{\left({\left(f_h^t\right)}_j\right)}^2+{\left({\left(f_v^t\right)}_j\right)}^2}},$ j＝1，…，N

14) the middle filtering image that this is obtained Assignment is given intermediate images

, and the value of t increased by 1,

15) if the value of t more than or equal to filtering total degree t _Max, then carried out for the 16th step, otherwise returned for the 11st step,

16) with intermediate images

Save as filtering image

, with this filtering image

As initial pictures, the corresponding increase by 1 of the value of m turned back to for the 7th step, repeated this process and reached the total degree I of iteration up to m, can obtain the image of final imaging.

Images