CN103116874B - A kind of picture breakdown method and system of maximum extracted noise useful information - Google Patents

Abstract

The invention discloses a kind of picture breakdown method, be specially: (1) initial decomposition original image I obtains initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾; (2) with initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V; (3) final smooth image U and final noise image V is added acquisition restored image G.The present invention also provides the system realizing above-mentioned decomposition method.Smooth part and noise section in Simultaneous Extracting Image of the present invention, and the useful information in noise section is incorporated into smooth image, obtain more genuine and believable restored image, thus make follow-up image procossing as more accurate in feature extraction, Iamge Segmentation and tissue typing etc.

Description

A kind of picture breakdown method and system of maximum extracted noise useful information

Technical field

The invention belongs to Digital Image Processing and medical imaging technology crossing domain, particularly a kind of picture breakdown method and system of maximum extracted noise useful information.

Background technology

Along with computer technology to be combined more and more closely with other kinds field subject and common development in all fields, no matter the technology and the means that are applied to medical image acquisition and analysis all have considerable raising from kind quantity or from practical function, this also make more directly perceived compared to TCM diagnosis method, accurately, the application of the higher medical image auxiliary diagnosis of efficiency in clinical be more extensive, therefore clinical diagnosis also becomes more accurately, convenient.

The development of medical image technology can be divided into the innovation of image acquisition hardware and improvement and the improvement to later image disposal route.Owing to (comprising the ultrasonic of comparatively early invention based on current any medical image acquisition means, computed tomography Computed tomography, CT) technology is to Magnetic resonance imaging (Magnatic ResonanceImaging) in recent years, positron emission tomography image (PositronEmission Tomography, PET) image etc.) obtained all can not be accomplished to can't harm completely (wherein more or less there is picture noise), therefore image reconstruction clinical diagnosis being had to material impact is being carried out, feature extraction, before the steps such as rim detection, very necessary to the image pre-service carried out to a certain degree.In traditional image procossing, being usually used as by the noise in the image obtained at first is useless interfere information, by adopting various filtering algorithm, as the denoising methods such as total variation denoising are removed.But research in recent years finds, the noise of medical image also contains useful information to a certain extent, as: at ultrasound field, Michailovich, O.and Tannenbaum, A., " De-speckling ofultrasound images, " IEEE Trans.Ultrason., Ferroelec., Freq.Contr.53 (1), 64 – 78 (2006).Therefore, before image post-processed, image is decomposed, instead of simple denoising, be necessary.And the theoretical background of picture breakdown has been proved; as: Meyer; Y.; [Oscillating patterns in image processing and nonlinear evolution equations]; vol.22 ofUniversity Lecture Series; American Mathematical Society, Providence, RI (2001); Osher, S., Sole, A., and Vese, L., " Image decomposition and restoration using total variationminimization and the H-1norm, " Multiscale Modeling and Simulation1,349370 (2003) etc., the use Bayesian frame maximum a posteriori probability (MAP) that the present invention proposes is estimated, the picture breakdown method for maximum extracted useful information from noise is a kind of selection of worth trial.

Summary of the invention

The speckle noise part existed is removed for the simple denoising method of employing of the prior art, likely lose the problems such as the useful information wherein comprised, the invention provides a kind of picture breakdown method and system of maximum extracted noise useful information, smooth part and noise section in Simultaneous Extracting Image, useful information in noise section is incorporated into smooth image, obtains more genuine and believable restored image.

A kind of picture breakdown method, is specially:

(1) initial decomposition original image I obtains initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

(2) with initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V;

(3) final smooth image U and final noise image V is added acquisition restored image G.

Further, described step (2) adopts alternating minimization method or least square method.

Further, described step (2) is specially:

(201) initialization iterations k=0, extracts initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾;

(202) smooth image of kth+1 time is calculated $U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{{\▿U}^{(k+1)}}{\left|\right|{\▿U}^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, represent gradient, || || represent norm;

(203) noise image of kth+1 time is calculated

V ^(k+1)＝V ^(k)+βp ^(k), $p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

${\▿}_V{E_V}^{\left(k\right)}=2(U^{\left(k\right)}+V^{\left(k\right)}-I)+{\mathrm{\λ}}_2[2\exp \left({2V}^{\left(k\right)}\right)-2],$ ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left({2V}^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

(204) if meet simultaneously $\frac{\underset{\mathrm{\Ω}}{\∫}(U^{(k+1)}-U^{\left(k\right)})\mathrm{dx}}{\underset{\mathrm{\Ω}}{\∫}{U}^{\left(k\right)}\mathrm{dx}}\≤U_{\mathrm{\ϵ}}$ With $\frac{\underset{\mathrm{\Ω}}{\∫}(V^{(k+1)}-V^{\left(k\right)})\mathrm{dy}}{\underset{\mathrm{\Ω}}{\∫}{V}^{\left(k\right)}\mathrm{dy}}\≤V_{\mathrm{\ϵ}},$ Then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise enter step (205), Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value;

(205) if iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise upgrade iterations k=k+1, return step (202).

Further, described step (2) is specially:

(211) initialization iterations k=0, extracts initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾;

(212) noise image of kth+1 time is calculated

V ^(k+1)＝V ^(k)+βp ^(k), $p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

${\▿}_V{E_V}^{\left(k\right)}=2(U^{\left(k\right)}+V^{\left(k\right)}-I)+{\mathrm{\λ}}_2[2\exp \left({2V}^{\left(k\right)}\right)-2],$ ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left({2V}^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

(213) smooth image of kth+1 time is calculated $U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{{\▿U}^{(k+1)}}{\left|\right|{\▿U}^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, represent gradient, || || represent norm;

(214) if meet simultaneously $\frac{\underset{\mathrm{\Ω}}{\∫}(U^{(k+1)}-U^{\left(k\right)})\mathrm{dx}}{\underset{\mathrm{\Ω}}{\∫}{U}^{\left(k\right)}\mathrm{dx}}\≤U_{\mathrm{\ϵ}}$ With $\frac{\underset{\mathrm{\Ω}}{\∫}(V^{(k+1)}-V^{\left(k\right)})\mathrm{dy}}{\underset{\mathrm{\Ω}}{\∫}{V}^{\left(k\right)}\mathrm{dy}}\≤V_{\mathrm{\ϵ}},$ Then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise enter step (215), Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value;

(215) if iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise upgrade iterations k=k+1, return step (212).

Further, described step (1) adopts any one in experience setting method, simple initialization method, random value initial method and total variation Denoising Algorithm.

A kind of picture breakdown system, comprises

First module, obtains initial smooth image U for initial decomposition original image I ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

Second module, for initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V;

3rd module, obtains restored image G for being added by final smooth image U and final noise image V.

Further, described second module adopts alternating minimization method or least square method.

Further, described second module comprises:

201st submodule, for initialization iterations k=0, extracts initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾;

202nd submodule, for calculating the smooth image of kth+1 time

$U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{{\▿U}^{(k+1)}}{\left|\right|{\▿U}^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, represent gradient, || || represent norm;

203rd submodule, for calculating the noise image of kth+1 time

V ^(k+1)＝V ^(k)+βp ^(k), $p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

${\▿}_V{E_V}^{\left(k\right)}=2(U^{\left(k\right)}+V^{\left(k\right)}-I)+{\mathrm{\λ}}_2[2\exp \left({2V}^{\left(k\right)}\right)-2],$ ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left({2V}^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

204th submodule, if meet for judging simultaneously with then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter the 3rd module, otherwise enter the 205th submodule, Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value;

205th submodule, if for judging that iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter the 3rd module, otherwise upgrade iterations k=k+1, return the 202nd submodule.

Further, described second module comprises:

211st submodule, for initialization iterations k=0, extracts initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾;

212nd submodule, for calculating the noise image of kth+1 time

V ^(k+1)＝V ^(k)+βp ^(k), $p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

${\▿}_V{E_V}^{\left(k\right)}=2(U^{\left(k\right)}+V^{\left(k\right)}-I)+{\mathrm{\λ}}_2[2\exp \left({2V}^{\left(k\right)}\right)-2],$ ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left({2V}^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

213rd submodule, for calculating the smooth image of kth+1 time $U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{{\▿U}^{(k+1)}}{\left|\right|{\▿U}^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, represent gradient, || || represent norm;

214th submodule, if meet for judging simultaneously with then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise enter the 215th submodule, Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value;

215th submodule, if for judging that iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter the 3rd module, otherwise upgrade iterations k=k+1, return the 212nd submodule.

Further, any one in described first module employing experience setting method, simple initialization method, random value initial method and total variation Denoising Algorithm carries out initial decomposition.

In sum, the present invention is directed to the existing disposal route Problems existing to noise in image spot, propose a kind of picture breakdown method and system of maximum extracted noise useful information, it takes full advantage of the noise section useful information of image, be integrated in the exploded view picture finally obtained, instead of simply carry out filtering removal speckle noise to image, thus make follow-up image procossing as more accurate in feature extraction, Iamge Segmentation and tissue typing etc.

Further, take full advantage of the priori of image, by signal and noise model, the iterative step of a kind of maximum extracted noise useful information derived under maximum a posteriori probability framework, it makes use of the method for the minimum optimization functional of iteration, simplifies calculating.

Accompanying drawing explanation

Fig. 1 is the picture breakdown method flow diagram of maximum extracted noise useful information of the present invention;

Fig. 2 is the result figure of embodiment one; Wherein, Fig. 2 (a) is original image; Fig. 2 (b) is initialization smooth image; Fig. 2 (c) is initialization noise image; Fig. 2 (d) is final smooth image; Fig. 2 (e) is final noise image; Fig. 2 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 2 (g) is residual image;

Fig. 3 is the sub-similarity test result figure of the inside and outside Bart Charlie of " grey value characteristics " arteria carotis implementing sample one; Wherein, Fig. 3 (a) is original image gray scale characteristic image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 3 (b) is the inside and outside experience density of original image gray scale arteria carotis and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis); Fig. 3 (c) is smooth image grey value characteristics image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 3 (d) is experience density inside and outside smooth image gray-scale value arteria carotis and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis);

Fig. 4 is the sub-similarity test result figure of the inside and outside Bart Charlie of " local entropy feature " arteria carotis implementing sample one: wherein, Fig. 4 (a) is original image local entropy characteristic image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 4 (b) is the inside and outside experience density of original image local entropy arteria carotis and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis); Fig. 4 (c) is smooth image local entropy characteristic image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 4 (d) is experience density inside and outside smooth image local entropy arteria carotis and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis);

Fig. 5 is the sub-similarity test result of the inside and outside Bart Charlie of " in local value tag " arteria carotis implementing sample one: wherein, Fig. 5 (a) is original image local intermediate value characteristic image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 5 (b) is the original image local inside and outside experience density of intermediate value arteria carotis and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis); Fig. 5 (c) is smooth image local intermediate value characteristic image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 5 (d) is experience density inside and outside the intermediate value arteria carotis of smooth image local and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis);

Fig. 6 is that embodiment one changes initialization mode " random initializtion (U ⁽⁰⁾, V ⁽⁰⁾random one can not equal 0 value) " result figure: wherein, Fig. 6 (a) is original image; Fig. 6 (b) is initialization smooth image; Fig. 6 (c) is initialization noise image; Fig. 6 (d) is final smooth image; Fig. 6 (e) is final noise image; Fig. 6 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 6 (g) is residual image;

Fig. 7 is that embodiment changes variable element setting " parameter lambda ₁and λ ₂equal respectively to be 2 and 1, iterations equals 9 " result figure: wherein, Fig. 7 (a) is original image; Fig. 7 (b) is initialization smooth image; Fig. 7 (c) is initialization noise image; Fig. 7 (d) is final smooth image; Fig. 7 (e) is final noise image; Fig. 7 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 7 (g) is residual image;

Fig. 8 is that embodiment changes variable element setting " parameter lambda ₁and λ ₂equal respectively to be 5 and 5, iterations equals 5 " result figure: wherein, Fig. 8 (a) is original image; Fig. 8 (b) is initialization smooth image; Fig. 8 (c) is initialization noise image; Fig. 8 (d) is final smooth image; Fig. 8 (e) is final noise image; Fig. 8 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 8 (g) is residual image;

Fig. 9 is result figure: Fig. 9 (a) of enforcement sample two is original image; Fig. 9 (b) is initialization smooth image; Wherein, Fig. 9 (c) is initialization noise image; Fig. 9 (d) is final smooth image; Fig. 9 (e) is final noise image; Fig. 9 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 9 (g) is residual image;

Figure 10 is result figure: Figure 10 (a) of enforcement sample three is original image; Figure 10 (b) is initialization smooth image; Wherein, Figure 10 (c) is initialization noise image; Figure 10 (d) is final smooth image; Figure 10 (e) is final noise image; Figure 10 (f) is restored image (superposing final smooth image and final resolution noise image); Figure 10 (g) is residual image;

Figure 11 is result figure: Figure 11 (a) of enforcement sample four is original image; Figure 11 (b) is initialization smooth image; Wherein, Figure 11 (c) is initialization noise image; Figure 11 (d) is final smooth image; Figure 11 (e) is final noise image; Figure 11 (f) is restored image (superposing final smooth image and final resolution noise image); Figure 11 (g) is residual image.

Figure 12 is result figure: Figure 12 (a) of enforcement sample five is original image; Figure 12 (b) is initialization smooth image; Wherein, Figure 12 (c) is initialization noise image; Figure 12 (d) is final smooth image; Figure 12 (e) is final noise image; Figure 12 (f) is restored image (superposing final smooth image and final resolution noise image); Figure 12 (g) is residual image.

Embodiment

In order to make object of the present invention, technical scheme and advantage clearly understand, below in conjunction with drawings and Examples, the present invention is further elaborated.Should be appreciated that specific embodiment described herein only in order to explain the present invention, be not intended to limit the present invention.

Succinct in order to compose a piece of writing, arrange following noun and be called for short:

Original image I

Smooth image U

Noise image V

Iteration smooth image initial value U ⁽⁰⁾

Iteration noise image initial value V ⁽⁰⁾

The smooth image U that kth time iteration obtains ^(k)

The noise image V that kth time iteration obtains ^(k)

Restored image G=U+V

Residual image W=I-(U+V)

The minimum optimization problem of energy functional represents argnin

Technical thought of the present invention is: largest optimization problem picture breakdown being considered as a posterior probability, namely bayesian criterion and negative logarithm operation is utilized, be converted into E (U, the minimum optimization problem of V)=-logP (I|U, V)-logP (U)-logP (V).In E (U, V), Section 1 describes residual image W=I-(U+V), and Section 2 describes smooth image U, and Section 3 describes noise image V.Thus, can obtain the minimum optimization problem of energy functional, solve this variational problem and obtain smooth image U and noise image V, complete picture breakdown task.Finally, the noise image V after the smooth image U Sum decomposition after decomposition is added, obtains the restored image G of original image.Noise image V after smooth image U Sum decomposition after decomposition, can carry out feature extraction respectively, obtain respective characteristic information, and then for follow-up image procossing (as Iamge Segmentation, classification, identification etc.).

Picture breakdown process flow diagram of the present invention, as shown in Figure 1.The derivation of core procedure of the present invention as shown in Figure 2.The picture breakdown method of a kind of maximum extracted noise useful information provided by the invention, comprises the following steps:

(1) initial decomposition original image I obtains initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾.

This step can adopt experience setting method (according to the given initial value of user's experience), simple initial method (directly makes U ⁽⁰⁾=I, V ⁽⁰⁾=I), random value initial method (U ⁽⁰⁾, V ⁽⁰⁾random one can not equal 0 value), total variation denoising model method, the present invention illustrates for the total variation denoising model of classics.

Total variation denoising model (also known as Rudin-Osher-Fatemi(ROF) model) from original image I, recover smooth image U by minimization of energy functional.

$\underset{U,V=I-U}{\arg \min }{\∫}_{\mathrm{\Ω}}|\▿U|+\mathrm{\λ}{\left|\right|V\left|\right|}_2^2---\left(1\right)$

Wherein, Ω is plane R ²bounded open subset, and V=I-U.Therefore (1) formula also can be expressed as:

$\underset{U}{\arg \min }{\∫}_{\mathrm{\Ω}}|\▿U|+\mathrm{\λ}{\left|\right|I-U\left|\right|}_2^2---\left(2\right)$

Smooth image U belongs to bounded variation (Bounded Variation, BV) function space, noise image V belongs to L2 norm space (on the function space that the quadractically integrable function on this measure space is formed, then can define L2 norm, norm is defined as: the square root of the integrated square of the absolute value of function.In addition L2 space can also define inner product.In other words, L2 (E) is exactly the entirety of the function that integrated square is limited on E).

In energy functional formula (1), Section 1 is regularization term, is the total variance norm (TV energy) of image, depends on the variation amplitude of image; Section 2 is item true to nature, which control the difference of smooth image U and original image I, and representing that smooth image U is to the matching degree of original image I, is the L of texture/noise V ²norm.λ > 0, is generally taken at [0,1] interval.λ is Lagrange (Lagrange) multiplier weighting coefficient, and its role is to keep the balance in energy functional formula between regularization term and item true to nature, the change of λ value can cause energy functional minimal value to obtain different results.Total variation model is a convex function, so its solution has existence and uniqueness.

Formula (1) also can be expressed as: $\underset{f=u+v,(u,v)\∈{\mathrm{BV}\×L}^2}{\inf }{\∫}_{\mathrm{\Ω}}|\▿u|+\mathrm{\λ}{\left|\right|f-u\left|\right|}_{L^2\left(\mathrm{\Ω}\right)}^2$

The great advantage of total variation denoising ROF model is, it while removal noise, can retain the edge of image.

In the present invention, the present invention defines the initial value U of U ⁽⁰⁾for total variation denoising Rudin-Osher-Fatemi(ROF) minimum value of model, that is:

$U^{\left(0\right)}=\underset{U}{\arg \min }{\∫}_{\mathrm{\Ω}}|\▿U|+\mathrm{\λ}{\left|\right|I-U\left|\right|}_2^2---\left(3\right)$

Wherein, λ is a scale parameter, is set to λ=0.3.The initial value of V is simply defined as original image I and U ⁽⁰⁾difference, that is: V ⁽⁰⁾=I-U ⁽⁰⁾.

(2) with initial noisc image volume and initial smooth image for initial value, be target to the maximum with the energy sum of smooth image U and noise image V, the final smooth image of iterative and noise image.

On the one hand after noise image change, smooth image also can respective change; On the other hand, image useful information is more, and energy is then larger, therefore maximizes extraction noise useful information target and can be converted into smooth image and noise image energy maximum target.This step can adopt the methods such as alternating minimization method, least square method to realize.The invention provides a kind of maximum extracted noise useful information step derived under maximum a posteriori probability framework, it makes use of the method for the minimum optimization functional of iteration, simplify calculating, the following detailed description of:

1) keep noise image V constant, the minimum optimization energy functional E of iteration _u.

Make V be constant, then have:

$E_u\left(U\right|V)={\left|\right|U+(I-V)\left|\right|}_2^2+{\mathrm{\λ}}_1{\left|\right|U\left|\right|}_{\mathrm{TV}}---\left(7\right)$

E described by formula (7) _u(U|V) function can be considered another ROF problem, so also can solve this problem under variation framework.

E _u(U|V) first variation as shown in the formula:

$\frac{{\mathrm{\δE}}_u\left(U\right|V)}{\mathrm{\δU}}=2((U+V)-I)-{\mathrm{\λ}}_1\mathrm{di}\left(v\frac{\▿U}{\left|\right|\▿U\left|\right|}\right)---\left(8\right)$

Adopt gradient descent algorithm, solved by the continuous iteration of following formula:

$\frac{\∂U}{\∂t}=-\frac{{\mathrm{\δE}}_u\left(U\right|V)}{\mathrm{\δU}}={\mathrm{\λ}}_1\mathrm{div}\left(\frac{\▿U}{\left|\right|\▿U\left|\right|}\right)+2(I-(U+V))---\left(9\right)$

Adopt Implicit numerical scheme, then in kth+1 iteration, U ^(k+1)change is shown below:

$U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{{\▿U}^{(k+1)}}{\left|\right|{\▿U}^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)}---\left(10\right)$

Wherein, div represents and asks divergence, represent gradient, || || represent norm, time parameter t can be a fixing empirical value or with the linear changing value of iterations.λ ₁be the weight coefficient of smooth part energy function, the general value 0 ~ 10 of weight coefficient, is adjusted according to the result run by user in practical application.

2) keep smooth part U constant, the minimum optimization energy functional E of iteration _v.

Make U be constant, corresponding modification done to E (U, V), then has:

E _v(V|U) function is convex function, and bis-times can be micro-about V.

This function about the first order derivative of V is:

${\▿}_V{E}_V=2(U+V-I)+{\mathrm{\λ}}_2[2\exp \left(2V\right)-2]---\left(12\right)$

The gloomy matrix in sea of second derivative is:

${\▿}_V^2{E}_V=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left(2V\right)\right)---\left(13\right)$

Wherein λ ₂for the weight coefficient of noise section energy function, the general value 0 ~ 10 of weight coefficient, is adjusted according to the result run by user in practical application.

The parameter needing user to adjust as required in the present invention has two, is the λ in functional formula E (U, V) ₁and λ ₂.

The gloomy matrix in sea of formula (13) is a diagonal matrix, therefore its inverse matrix can directly calculate.Just because of this, the present invention can solve E by Newton method _v(V|U) minimum optimization problem.

In kth+1 iteration, V alternates and is shown below:

V ^(k+1)＝V ^(k)+βp ^(k), $p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)}---\left(14\right)$

Wherein, ${\▿}_V{E_V}^{\left(k\right)}=2(U^{\left(k\right)}+V^{\left(k\right)}-I)+{\mathrm{\λ}}_2[2\exp \left({2V}^{\left(k\right)}\right)-2],$ ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left({2V}^{\left(k\right)}\right)\right),$ β represents along Newton direction p ^(k)carry out the step-length (empirical value, recommendation is 0.5 ~ 1) of linear search, exp represents index, and diag represents diagonal matrix.To constantly carry out such as formula the change described by (10) and formula (14), until when the relative change of U with V is very little, program is iteration convergence.

Wherein judge whether that the condition stopped is:

Smoothed image satisfies condition 1: wherein x is image U ^(k)pixel, Ω is real number space R ²bounded open subset, U _εbe the accuracy requirement be required to meet, accuracy requirement is higher, then U _εless

Noise image satisfies condition 2: wherein y is image V ^(k)pixel, Ω is real number space R ²bounded open subset, V _εbe the accuracy requirement be required to meet, accuracy requirement is higher, then V _εless.

Iterated conditional satisfies condition 3: iterations is less than or equal to 100;

The iteration order of U and V, little to Influence on test result.Constant with first U below, then the iterative process of this step of the constant detailed description of V:

(21) initialization iterations k=0, extracts initial noisc image U ⁽⁰⁾with initial smooth image V ⁽⁰⁾;

(22) by smooth image U that kth time iteration obtains ^(k)with the noise image V that kth time iteration obtains ^(k)substitution formula (10), obtains the smooth image after kth+1 iteration;

(23) noise image iteration: the smooth image U that kth+1 iteration is obtained ^(k+1)with the noise image V that kth time iteration obtains ^(k)substitution formula (14), obtains the noise image V after kth+1 iteration ^(k+1).

(24) judge whether to stop: compare (x is image U ^(k)pixel) with (y is image V ^(k)pixel).If the two is all set up, then iteration ends, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3); Otherwise, enter step (D).

(25) if iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1)iteration, enters step (3), otherwise upgrades iterations k=k+1, returns step (22).

Iterations of the present invention changes according to the difference of image type.General, iterations is no more than 100, otherwise just needs to carry out parameter adjustment to adapt to image request.The smooth image that decomposition obtains and noise image can be used for successive image analysis, carry out the process such as rim detection and feature extraction with support study dies personnel to it.

(3) smooth image and noise image addition acquisition restored image that obtain will finally be decomposed.

Being added finally decomposing the smooth image U that obtains and noise image V, obtaining restored image G, i.e. G=U+V; Deduct smooth image U and noise image V with original image I, namely obtain residual image W, be i.e. W=I-(U+V).

Enforcement of the present invention can be applicable to ultrasonoscopy, visible images, Magnetic resonance imaging (Magnatic ResonanceImaging), confocal images (Confocal microscope image), computed tomography images Computed tomography, CT), CT angiography image (CT angiography, CTA), ultrasound computed tomography image (Ultrasound CT, UCT), fluorescence labeling image (Fluorescent marker image), positron emission tomography image (Positron Emission Tomography, PET) etc., various dissimilar image.

By reference to the accompanying drawings, below, the present invention introduces its concrete implementation method for five kinds of dissimilar images (ultrasonoscopy, nuclear-magnetism image, fluorescence labeling image, confocal images and MIcrosope image) by reference to the accompanying drawings.

Embodiment one ultrasonoscopy decomposes (Fig. 2, Fig. 6, Fig. 7, Fig. 8)

Embodiment of the present invention dividing method, comprises following four steps:

(1) be loaded into original image Fig. 2 (a) to be decomposed, carry out total variation initialization, parameter lambda=0.3, obtain initial smooth image Fig. 2 (b) and initial noisc image graph 2(c).Noise in noise image obeys Fisher-Tippett distribution.

(2) adopt the minimum optimization energy Functional Approach of alternating iteration, work as parameter lambda ₁and λ ₂when equaling respectively to be 5 and 3, and when iterations equals 10, meet the condition of convergence.U in this example _ε=0.1%, V _ε=0.1%.After iteration completes, the final smooth image of decomposition is as shown in Fig. 2 (d), and the final noise image of decomposition is as shown in (e) in Fig. 2.

(3) will finally decompose the smooth image and noise image addition acquisition restored image that obtain, restored image is as shown in Fig. 2 (f), and residual image is as shown in (g) image in Fig. 2.

The smooth image that decomposition obtains and noise image can be used for the analysis of successive image feature extraction, carry out the process such as rim detection and feature extraction with support study dies personnel to it.Feature (referring to Fig. 3, Fig. 4 and Fig. 5) such as " gray-scale value ", " local entropy " and " local intermediate value " of relatively this example original image and final smooth image below:

The ultrasonic carotid images of the embodiment of the present invention one Fig. 2 (a) is example, introduces Bhattacharyya similarity concept, judges the resolution characteristic between artery and background by the Bhattacharyya similarity of experience density.If the Bhattacharyya distance between artery experience density and background experience density is larger, key diagram is better as the separating capacity between artery and background.

Here, in order to easy, directly adopt the gray probability intensity profile empirically density in region.Specific practice is: the characteristic image first asking for certain feature, then asks for the gray probability distribution of background area and target area in this characteristic image, recycling formula calculate the sub-similarity of Bart Charlie between these two intensity profile.Utilize the sub-similarity of this Bart Charlie can judge the ability in this feature differentiation target and background region.

Here, the present invention compares the features such as " gray-scale value ", " local entropy " and " local intermediate value " of original image and final smooth image, then calculates and compare the sub-similarity of Bart Charlie of each feature experience density in artery and outside artery respectively.

Fig. 3 is the sub-similarity test result of " grey value characteristics " arteria carotis inside and outside Bart Charlie of the original image (Fig. 2 (a)) of embodiment one: Fig. 3 (a), Fig. 3 (c) are respectively original image, the gray-value image of smooth image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 3 (b), Fig. 3 (d) are respectively the sub-similarity of original image, the inside and outside experience density of smooth image gray-scale value arteria carotis and Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis);

Fig. 4 be the original image (Fig. 2 (a)) of embodiment one " local entropy feature " arteria carotis inside and outside the sub-similarity test result of Bart Charlie: Fig. 4 (a), Fig. 4 (c) are respectively original image, the local entropy image of smooth image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 4 (b), Fig. 4 (d) are respectively the sub-similarity of original image, the inside and outside experience density of local entropy arteria carotis of smooth image and Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis);

Fig. 5 be the original image (Fig. 2 (a)) of embodiment one " local in value tag " arteria carotis inside and outside the sub-similarity test result of Bart Charlie; Fig. 5 (a), Fig. 5 (c) are respectively original image, the local median image of smooth image and interior epicardium contours (black dotted lines represents arteria carotis epicardium contours, and solid white line represents carotid artery intima profile); Fig. 5 (b), Fig. 5 (d) are respectively the original image local inside and outside experience density of intermediate value arteria carotis and the sub-similarity of Bart Charlie (black color dots line represents the inner experience density of arteria carotis, and grey filled lines represents intermembranous endless belt experience density inside and outside arteria carotis);

As can be seen from Fig. 3, Fig. 4, Fig. 5, inside and outside the gray-scale value of smooth image, local entropy and local intermediate value arteria carotis, Bhattacharyya similarity is all lower than the original image of its correspondence.This illustrates, the smooth image after decomposition is compared to original image, and its gray-scale value, local entropy and local intermediate value feature are more conducive to distinguishing arteria carotis and background.

Embodiment one also only changes " initialization " mode on the basis of above-mentioned steps Fig. 2 (a), and " total variation initialization " is replaced with " random initializtion (U ⁽⁰⁾, V ⁽⁰⁾random one can not equal 0 value) " (Fig. 6), and other steps are constant.Wherein: Fig. 6 (a) is original image; Fig. 6 (b) is initialization smooth image; Fig. 6 (c) is initialization noise image; Fig. 6 (d) is final smooth image; Fig. 6 (e) is final noise image; Fig. 6 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 6 (g) is residual image.In addition, also with " simple initialization (directly makes U ⁽⁰⁾=I, V ⁽⁰⁾=I) " compare.Experiment proves: random initializtion and simple initialization also can be decomposed and obtained more real restored image, but total variation denoising initialization effect is better than simple initialization and random initializtion method.

Embodiment one also only changes " parameter " and arranges, by " parameter lambda on the basis of above-mentioned steps ₁and λ ₂equal respectively to be 5 and 3, iterations equals 10 " replace with " parameter lambda ₁and λ ₂equal respectively to be 2 and 1, iterations equals 9 " (Fig. 7), " parameter lambda ₁and λ ₂equal respectively to be 5 and 5, iterations equals 5 " (Fig. 8), and other steps are constant.Wherein: parameter lambda ₁and λ ₂equaling respectively to be 2 and 1, when iterations equals 9, Fig. 7 (a) is original image; Fig. 7 (b) is initialization smooth image; Fig. 7 (c) is initialization noise image; Fig. 7 (d) is final smooth image; Fig. 7 (e) is final noise image; Fig. 7 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 7 (g) is residual image; Parameter lambda ₁and λ ₂equaling respectively to be 5 and 5, when iterations equals 5, Fig. 8 (a) is original image; Fig. 8 (b) is initialization smooth image; Fig. 8 (c) is initialization noise image; Fig. 8 (d) is final smooth image; Fig. 8 (e) is final noise image; Fig. 8 (f) is restored image (superposing final smooth image and final resolution noise image); Fig. 8 (g) is residual image.

Experiment proves: after decomposing, the flatness of smooth image is fine, and remains edge and the details of original image preferably.Therefore Ultrasonic Image Denoising can be effectively applied to the picture breakdown algorithm that invention proposes, different to the demand of image denoising degree according to user, reasonably adjust parameter lambda ₁and λ ₂parameter value.In addition, suitably λ is turned down ₂parameter, time namely less demanding to the matching degree of noise section and its distribution, this invention also can be used for the denoising of other ordinary numbers image, and is not limited to ultrasonoscopy.In the present invention, iterations is generally within 100.

As can be seen from Fig. 6, Fig. 7, Fig. 8, different initial methods, different parameters are impartial has certain influence to decomposition final effect, and concrete setting should combine with object with application background.

Embodiment two nuclear-magnetism picture breakdown (Fig. 9)

(1) be loaded into original image Fig. 9 (a) to be decomposed, carry out total variation initialization, obtain initial smooth image Fig. 9 (b) and initial noisc image graph 9(c).Noise Rayleigh distributed in noise image.

(2) adopt least square method, work as parameter lambda ₁and λ ₂when equaling respectively to be 5 and 3, and when iterations equals 52, meet the condition of convergence.After iteration completes, the final smooth image of decomposition is as shown in Fig. 9 (d), and the final noise image of decomposition is as shown in (e) in Fig. 9.

(3) will finally decompose the smooth image and noise image addition acquisition restored image that obtain, restored image is as shown in Fig. 9 (f), and the residual image that the present invention obtains is as shown in (g) image in Fig. 9.

Embodiment three fluorescence labeling picture breakdown (Figure 10)

(1) be loaded into original image Figure 10 (a) to be decomposed, carry out random initializtion, obtain initial smooth image Figure 10 (b) and initial noisc image graph 10(c).Noise in noise image obeys Na Kajiami (Nakagami) distribution.

(2) adopt the minimum optimization energy Functional Approach of alternating iteration, work as parameter lambda ₁and λ ₂when equaling respectively to be 5 and 3, and when iterations equals 15, meet the condition of convergence.After iteration completes, the final smooth image of decomposition is as shown in Figure 10 (d), and the final noise image of decomposition is as shown in (e) in Figure 10.

(3) will finally decompose the smooth image and noise image addition acquisition restored image that obtain, restored image is as shown in Figure 10 (f), and residual image is as shown in (g) image in Figure 10.

Embodiment four confocal images decomposes (Figure 11)

(1) be loaded into original image Figure 11 (a) to be decomposed, carry out simple initialization, obtain initial smooth image Figure 11 (b) and initial noisc image graph 11(c).Noise in noise image obeys the distribution of card side.

(2) adopt the minimum optimization energy Functional Approach of alternating iteration, work as parameter lambda ₁and λ ₂when equaling respectively to be 5 and 3, and when iterations equals 15, meet the condition of convergence.After iteration completes, the final smooth image of decomposition is as shown in Figure 11 (d), and the final noise image of decomposition is as shown in (e) in Figure 11.

(3) will finally decompose the smooth image and noise image addition acquisition restored image that obtain, restored image is as shown in Figure 11 (f), and residual image is as shown in (g) image in Figure 11.

Embodiment five MIcrosope image decomposes (Figure 12)

(1) be loaded into original image Figure 12 (a) to be decomposed, carry out total variation initialization, obtain initial smooth image Figure 12 (b) and initial noisc image graph 12(c).Noise Gaussian distributed in noise image.

(2) adopt the minimum optimization energy Functional Approach of alternating iteration, work as parameter lambda ₁and λ ₂when equaling respectively to be 5 and 3, and when iterations equals 15, meet the condition of convergence.After iteration completes, the final smooth image of decomposition is as shown in Figure 12 (d), and the final noise image of decomposition is as shown in (e) in Figure 12.

(3) will finally decompose the smooth image and noise image addition acquisition restored image that obtain, restored image is as shown in Figure 12 (f), and residual image is as shown in (g) image in Fig. 2.

Those skilled in the art will readily understand; the foregoing is only preferred embodiment of the present invention; not in order to limit the present invention, all any amendments done within the spirit and principles in the present invention, equivalent replacement and improvement etc., all should be included within protection scope of the present invention.

Claims (4)

1. a picture breakdown method, is characterized in that, is specially:

(1) initial decomposition original image I obtains initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

(2) with initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V;

(3) final smooth image U and final noise image V is added acquisition restored image G;

Described step (2) is specially:

(201) initialization iterations k=0, extracts initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

(202) smooth image of kth+1 time is calculated $U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{\▿U^{(k+1)}}{\left|\right|\▿U^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, and ▽ represents gradient, || || represent norm;

(203) noise image of kth+1 time is calculated

$V^{(k+1)}=V^{\left(k\right)}+\mathrm{\β}{p}^{\left(k\right)},p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

▽ _VE _V ^(k)＝2(U ^(k)+V ^(k)-I)+λ ₂[2exp(2V ^(k))-2]， ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left(2V^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

(204) if meet simultaneously $\frac{\underset{\mathrm{\Ω}}{\∫}(U^{(k+1)}-U^{\left(k\right)})\mathrm{dx}}{\underset{\mathrm{\Ω}}{\∫}{U}^{\left(k\right)}\mathrm{dx}}\≤U_{\mathrm{\ϵ}}$ With $\frac{\underset{\mathrm{\Ω}}{\∫}(V^{(k+1)}-V^{\left(k\right)})\mathrm{dy}}{\underset{\mathrm{\Ω}}{\∫}{V}^{\left(k\right)}\mathrm{dy}}\≤V_{\mathrm{\ϵ}},$ Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise, enter step (205);

(205) if iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise upgrade iterations k=k+1, return step (202).

2. a picture breakdown method, is characterized in that, is specially:

(1) initial decomposition original image I obtains initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

(2) with initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V;

(3) final smooth image U and final noise image V is added acquisition restored image G;

Described step (2) is specially:

(211) initialization iterations k=0, extracts initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

(212) noise image of kth+1 time is calculated

$V^{(k+1)}=V^{\left(k\right)}+\mathrm{\β}{p}^{\left(k\right)},p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

▽ _VE _V ^(k)＝2(U ^(k)+V ^(k)-I)+λ ₂[2exp(2V ^(k))-2]， ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left(2V^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

(213) smooth image of kth+1 time is calculated $U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{\▿U^{(k+1)}}{\left|\right|\▿U^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, and ▽ represents gradient, || || represent norm;

(214) if meet simultaneously $\frac{\underset{\mathrm{\Ω}}{\∫}(U^{(k+1)}-U^{\left(k\right)})\mathrm{dx}}{\underset{\mathrm{\Ω}}{\∫}{U}^{\left(k\right)}\mathrm{dx}}\≤U_{\mathrm{\ϵ}}$ With $\frac{\underset{\mathrm{\Ω}}{\∫}(V^{(k+1)}-V^{\left(k\right)})\mathrm{dy}}{\underset{\mathrm{\Ω}}{\∫}{V}^{\left(k\right)}\mathrm{dy}}\≤V_{\mathrm{\ϵ}},$ Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise, enter step (215);

(215) if iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise upgrade iterations k=k+1, return step (212).

3. a picture breakdown system, comprises

First module, obtains initial smooth image U for initial decomposition original image I ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

Second module, for initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V;

3rd module, obtains restored image G for being added by final smooth image U and final noise image V;

Described second module comprises:

201st submodule, for initialization iterations k=0, extracts initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

202nd submodule, for calculating the smooth image of kth+1 time

$U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{\▿U^{(k+1)}}{\left|\right|\▿U^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, and ▽ represents gradient, || || represent norm;

203rd submodule, for calculating the noise image of kth+1 time

$V^{(k+1)}=V^{\left(k\right)}+\mathrm{\β}{p}^{\left(k\right)},p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

▽ _VE _V ^(k)＝2(U ^(k)+V ^(k)-I)+λ ₂[2exp(2V ^(k))-2]， ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left(2V^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

204th submodule, if meet for judging simultaneously $\frac{\underset{\mathrm{\Ω}}{\∫}(U^{(k+1)}-U^{\left(k\right)})\mathrm{dx}}{\underset{\mathrm{\Ω}}{\∫}{U}^{\left(k\right)}\mathrm{dx}}\≤U_{\mathrm{\ϵ}}$ With $\frac{\underset{\mathrm{\Ω}}{\∫}(V^{(k+1)}-V^{\left(k\right)})\mathrm{dy}}{\underset{\mathrm{\Ω}}{\∫}{V}^{\left(k\right)}\mathrm{dy}}\≤V_{\mathrm{\ϵ}},$ Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter the 3rd module, otherwise, enter the 205th submodule;

205th submodule, if for judging that iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter the 3rd module, otherwise upgrade iterations k=k+1, return the 202nd submodule.

4. a picture breakdown system, is characterized in that, comprises

First module, obtains initial smooth image U for initial decomposition original image I ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

Second module, for initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾for initial value, be target to the maximum with the energy sum of smooth image and noise image, iterative final smooth image U and noise image V;

3rd module, obtains restored image G for being added by final smooth image U and final noise image V;

Described second module comprises:

211st submodule, for initialization iterations k=0, extracts initial smooth image U ⁽⁰⁾with initial noisc image V ⁽⁰⁾;

212nd submodule, for calculating the noise image of kth+1 time

$V^{(k+1)}=V^{\left(k\right)}+\mathrm{\β}{p}^{\left(k\right)},p^{\left(k\right)}=-{\left({\▿}_V^2{E_V}^{\left(k\right)}\right)}^{-1}\·{\▿}_V{E_V}^{\left(k\right)},$

▽ _VE _V ^(k)＝2(U ^(k)+V ^(k)-I)+λ ₂[2exp(2V ^(k))-2]， ${\▿}_V^2{E_V}^{\left(k\right)}=2I+4{\mathrm{\λ}}_2\mathrm{diag}\left(\exp \left(2V^{\left(k\right)}\right)\right),$

β represents along Newton direction p ^(k)carry out the step-length of linear search, exp represents index, and diag represents diagonal matrix;

213rd submodule, for calculating the smooth image of kth+1 time

$U^{(k+1)}=[{\mathrm{\λ}}_1\mathrm{div}\left(\frac{\▿U^{(k+1)}}{\left|\right|\▿U^{\left(k\right)}\left|\right|}\right)+2(I-(U^{\left(k\right)}+V^{\left(k\right)}))]t+U^{\left(k\right)},$ T is time parameter, λ ₁be the weight coefficient of smooth part energy function, div represents divergence, and ▽ represents gradient, || || represent norm;

214th submodule, if meet for judging simultaneously $\frac{\underset{\mathrm{\Ω}}{\∫}(U^{(k+1)}-U^{\left(k\right)})\mathrm{dx}}{\underset{\mathrm{\Ω}}{\∫}{U}^{\left(k\right)}\mathrm{dx}}\≤U_{\mathrm{\ϵ}}$ With $\frac{\underset{\mathrm{\Ω}}{\∫}(V^{(k+1)}-V^{\left(k\right)})\mathrm{dy}}{\underset{\mathrm{\Ω}}{\∫}{V}^{\left(k\right)}\mathrm{dy}}\≤V_{\mathrm{\ϵ}},$ Ω is plane real number space R ²bounded open subset, U _εand V _εfor accuracy requirement value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter step (3), otherwise, enter the 215th submodule;

215th submodule, if for judging that iterations k equals predetermined value, then final smooth image U=U ^(k+1)with noise image V=V ^(k+1), enter the 3rd module, otherwise upgrade iterations k=k+1, return the 212nd submodule.