CN103810731A - PET image reconstruction method based TV norm - Google Patents

Abstract

The invention discloses a PET (Positron Emission Tomography) image reconstruction method based a TV ( Total Variation) norm. The PET image reconstruction method is capable of reconstructing a PET image based on TV norm method by establishing a mathematical model of reconstruction problems and transforming an equivalence problem model, wherein in the process of reconstructing the PET image by virtue of a TV norm model, a TVAL3 algorithm is employed in the solving process of equivalence model problems. As a result, the TVAL3 algorithm is utilized effectively, and the problem of high requirement for data collection quantity in the reconstruction process of traditional PET image reconstruction is improved; experiments of comparing with the existing reconstruction method indicate that good reconstruction result can be maintained while the data collection quantity is reduced.

Description

A kind of PET image rebuilding method based on TV norm

Technical field

The invention belongs to PET technical field of imaging, be specifically related to a kind of based on TV(total variation) the PET image rebuilding method of norm.

Background technology

Positron emission tomography (Positron emission tomography, PET) be a kind of based on nuclear physics, molecular biological Medical Imaging Technology, it can be from molecular level the metabolic activity of observation of cell, for detection and the prevention of early stage disease provide effective foundation.First need in by the isotopically labeled infusion of medicine human body of radioactivity coordination carrying out PET when scanning, by blood circulation system, these materials will form certain distribution in each histoorgan in human body.Because the half life period of radioactivity coordination nucleic is shorter, and extremely unstable, to decay very soon, the positron producing in decay process and near free electron generation annihilation reaction, produce a pair of direction almost contrary, energy equate, size is the gamma photons of 511kev, via meet acquisition system to these paired photon with radiopharmaceutical distribution information process generate data for projection.By corresponding mathematical method, data for projection is carried out to inverting and solve, can reconstruct the spatial concentration distribution of the radiomaterial of human body.

In recent years, the increase of PET detector quantity, the raising of accuracy of data acquisition and clinically the demand of high resolving power, large area scope scan image has been proposed to new challenge to image reconstruction algorithm: the expansion of image dimension, to affect computing machine processing speed, even cause calculator memory deficiency, the advantage that in addition, must guarantee to make full use of image data amount aspect is to improve image resolution ratio.

At present, PET image rebuilding method is broadly divided into two classes: analytical method and iteration statistic law.Last class is mainly filtered back projection's method, and computing velocity is fast, but artifact is serious.For this reason, the iteration statistic law take maximum likelihood method as representative is suggested, because process of iteration is based on statistical models, good to fragmentary data adaptability, becomes gradually PET reconstruction algorithm research focus.Rebuild speed issue for solving, the extensive quick computing of matrix is constantly applied in this algorithm system, as order subset maximum likelihood method.And be directed to the statistical property of not considering prior estimate in algorithm, the maximum a posteriori estimation technique of being revised during posteriority is estimated has been proposed.In addition, introduced again least square method after statistical model is approximately to Gauss model, development is not with least square, the non-negative over-relaxation iterative method etc. of punishment weighting then.

State space system provides the new approaches of PET image reconstruction from a brand-new angle, by measure the uniform expression of equation and state equation according to actual Solve problems adjustment, to realize static state, dynamic reconstruction and prior estimate.Solve by correlation technique, as Kalman filtering, H ∞ filtering, particle filter etc. can adapt to different noisinesss, nonidentity operation condition, have certain advantage compared with traditional analytical method or iteration statistic law.But the state equation of its structure and measure in equation and relate to high-dimensional matrix stores and computing, and inversion process wherein easily causes low memory, is a large key factor of restriction state space algorithm feasibility.

Said method has higher requirement to the data volume gathering simultaneously, has also just caused current data acquisition time, and PET is longer sweep time, also requires the certain use amount of spike medicine.

Summary of the invention

For the existing above-mentioned technical matters of prior art, the invention provides a kind of PET image rebuilding method based on TV norm, solve the data volume that computing machine relates to gathering in the process of carrying out image reconstruction and had relatively high expectations, the problem that tracer agent consumption is large.

A PET image rebuilding method based on TV norm, comprises the steps:

(1) utilize detector to survey the biological tissue that is injected with radiomaterial, collect the coincidence counting vector of current time, and this coincidence counting vector is proofreaied and correct;

(2), according to PET image-forming principle, the measurement equation of setting up PET is as follows:

y=Dx

Wherein: D is system matrix, y is coincidence counting after proofreading and correct and for m dimensional vector, and x is PET CONCENTRATION DISTRIBUTION vector and is n-dimensional vector, and m and n are the natural number that is greater than 1;

(3) introduce TV norm by the measurement equation to described, the TV problem model that obtains PET is as follows:

$\begin{array}{cc}\underset{x}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left|\right|H_i\left(x\right)\left|\right|& s.t.y=\mathrm{Dx}\end{array}$

Wherein: H _i(x) for two element values in bivector and this bivector are respectively x _{i ,+1x}-x _iand x _{i ,+1y}-x _i; x _ifor the concentration value of i element in PET CONCENTRATION DISTRIBUTION vector x, x _{i ,+1x}for x in PET image corresponding to PET CONCENTRATION DISTRIBUTION vector x _ithe inner x that is expert at of institute _ithe concentration value of a rear column element, x _{i ,+1y}for x in PET image corresponding to PET CONCENTRATION DISTRIBUTION vector x _ix in column _ithe concentration value of a rear row element, i is natural number and 1≤i≤n;

(4) described TV problem model is solved, the enhancement mode Lagrangian function that obtains TV problem model is as follows:

$L_A(\mathrm{\ω},x)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\right|\left|{\mathrm{\ω}}_i\right||-v_i^T(H_i\left(x\right)-{\mathrm{\ω}}_i)+\frac{{\mathrm{\β}}_i}{2}{\left|\right|H_i\left(x\right)-{\mathrm{\ω}}_i\left|\right|}_2^2)-{\mathrm{\λ}}^T(\mathrm{Dx}-y)+\frac{\mathrm{\μ}}{2}{\left|\right|\mathrm{Dx}-y\left|\right|}_2^2$

Wherein: L _a(ω, x) is the enhancement mode Lagrangian function about ω and x, and ω is the gradient matrix corresponding with x, ω _ifor the two-dimensional gradient vector of corresponding i element in gradient matrix ω, v _iand β _ibe respectively the punishment vector sum penalty coefficient of corresponding i element, λ and μ are respectively the punishment vector sum penalty coefficient of Lagrangian function;

(5) described enhancement mode Lagrangian function is minimized and solved, estimate to obtain PET CONCENTRATION DISTRIBUTION vector x; And then carry out PET imaging according to the PET CONCENTRATION DISTRIBUTION vector x obtaining.

In described step (4), utilize TVAL3 algorithm (the TV restructing algorithm based on compressed sensing) to solve TV problem model, obtain the enhancement mode Lagrangian function of TV problem model.

In described step (5), by following iterative equation group, enhancement mode Lagrangian function is minimized and is solved:

x ^k+1=x ^k-α ^kZ ^k(x ^k)

$\begin{array}{cc}{\mathrm{\ω}}_i^{k+1}=\max \left\{\right||H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\β}}_i^k}||-\frac{1}{{\mathrm{\β}}_i^k},0\}*\frac{(H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\β}}_i^k})}{\left|\right|H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\β}}_i^k}\left|\right|}& {\mathrm{\β}}_i^k\≥{\mathrm{\β}}_i^{k-1}\end{array}$

Wherein: x ^kand x ^k+1be respectively the PET CONCENTRATION DISTRIBUTION vector of the k time iteration and the k+1 time iteration,

be the gradient matrix ω of the k+1 time iteration ^k+1the two-dimensional gradient vector of middle corresponding i element,

with

be respectively the penalty coefficient of corresponding i element in the k time iteration and the k-1 time iteration,

be the punishment vector of corresponding i element in the k time iteration, α ^kbe the middle transition vector of the k time iteration, Z ^k(x ^k) be about x ^kpenalty term function; H _i() is that two element values in bivector and this bivector are respectively h _{i ,+1x}-h _iand h _{i ,+1y}-h _i; Wherein, h _ifor i element value of vector in (), h _{i ,+1x}for h in the PET image that in (), vector is corresponding _ithe inner h that is expert at of institute _ia rear column element value, h _{i ,+1y}for h in the PET image that in (), vector is corresponding _ih in column _irear a line element value.

Described middle transition vector α ^kcalculation expression as follows:

${\mathrm{\α}}^k=\frac{{(x^k-x^{k-1})}^T(x^k-x^{k-1})}{{(x^k-x^{k-1})}^T(Z^k\left(x^k\right)-Z^k\left(x^{k-1}\right))}$

Wherein: Z ^k(x ^k-1) be about x ^k-1penalty term function, x ^k-1it is the PET CONCENTRATION DISTRIBUTION vector of the k-1 time iteration.

Described penalty term function Z ^k(x ^k) and Z ^k(x ^k-1) calculation expression as follows:

$Z^k\left(x^k\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[{\mathrm{\β}}_i^k{H}_i^T(-H_i\left(x^k\right)-{\mathrm{\ω}}_i^{k+1})-H_i^T\left(v_i^k\right)]+{\mathrm{\μ}}^k{D}^T({\mathrm{Dx}}^k-y)-D^T{\mathrm{\λ}}^k$

$Z^k\left(x^{k-1}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[{\mathrm{\β}}_i^k{H}_i^T(-H_i\left(x^{k-1}\right)-{\mathrm{\ω}}_i^k)-H_i^T\left(v_i^k\right)]+{\mathrm{\μ}}^k{D}^T({\mathrm{Dx}}^{k-1}-y)-D^T{\mathrm{\λ}}^k$

Wherein:

be the gradient matrix ω of the k time iteration ^kthe two-dimensional gradient vector of middle corresponding i element, λ ^kand μ ^kbe respectively the punishment vector sum penalty coefficient of Lagrangian function in the k time iteration;

for two element values in bivector and this bivector are respectively h _{i ,-1x}-h _iand h _{i ,-1y}-h _i; Wherein, h _ifor i element value of vector in (), h _{i ,-1x}for h in the PET image that in (), vector is corresponding _ithe inner h that is expert at of institute _iprevious column element value, h _{i ,-1y}for h in the PET image that in (), vector is corresponding _ih in column _iprevious row element value.

Described punishment vector

calculation expression as follows:

$v_i^k=v_i^{k-1}-{\mathrm{\β}}_i^{k-1}(H_i\left(x^{k-1}\right)-{\mathrm{\ω}}_i^{k-1})$

Wherein:

be the punishment vector of corresponding i element in the k-1 time iteration,

be the gradient matrix ω of the k-1 time iteration ^k-1the two-dimensional gradient vector of middle corresponding i element.

Described punishment vector λ ^kcalculation expression as follows:

λ ^k=λ ^k-1-μ ^k-1(Dx ^k-1-y) μ ^k≥μ ^k-1

Wherein: λ ^k-1and μ ^k-1be respectively the punishment vector sum penalty coefficient of Lagrangian function in the k-1 time iteration.

Carry out iterative computation by described iterative equation group, the PET CONCENTRATION DISTRIBUTION vector x obtaining as estimation when the PET CONCENTRATION DISTRIBUTION vector reaching after maximum iteration time or iteration convergence; Iteration convergence condition is as follows:

$\frac{{\left|\right|x^{k+1}-x^k\left|\right|}_2}{\left|\right|x^k\left|\right|}<\mathrm{\&rho;}$

Wherein: ρ is given convergence threshold, || || represent norm, || || ₂represent 2 norms.

The present invention effectively utilizes TVAL3 algorithm, has improved the problem that image data amount is had relatively high expectations existing in process of reconstruction based on conventional P ET image reconstruction; Relatively show with the experiment of existing method for reconstructing, reduce under prerequisite in image data amount, can keep rebuilding preferably effect.

Accompanying drawing explanation

Fig. 1 is the schematic flow sheet of PET image rebuilding method of the present invention.

Fig. 2 is the schematic flow sheet of alternating minimization module correction slack variable and PET image.

Fig. 3 (a) is the true value image about brain body mould.

Fig. 3 (b) is for adopting EM algorithm to rebuild the PET image of brain body mould.

Fig. 3 (c) is for adopting the inventive method to rebuild the PET image of brain body mould.

Fig. 4 is the inventive method and the EM algorithm comparison diagram about data reconstruction result.

Embodiment

In order more specifically to describe the present invention, below in conjunction with the drawings and the specific embodiments, technical scheme of the present invention is elaborated.

As shown in Figure 1, a kind of PET image rebuilding method based on TV norm, comprises step:

PET (positron emission tomography) scanner is surveyed the radiated signal sending in human body, through meeting and acquisition system processing, forms raw projections line, and deposits in hard disc of computer in the mode of sinogram; The sinogram that acquired original is arrived, take sinogram and known system matrix D as input item, calls correlation module.

S1. the principle of surveying according to PET is set up the basic model of Problems of Reconstruction;

S2. introduce TVAL3 separate with S1 in the optimization problem of model equivalence;

S3. initialization λ ⁰, μ ⁰; Initial point slack variable is set

and x ⁰;

S4. judge whether to meet iteration stopping condition, do not meet this condition and perform step S5, satisfied iteration stopping;

S5. arrange

with

from

with

start, calculate the slack variable in the Lagrange's equation that strengthens version according to alternating minimization method

with the image x rebuilding ^k+1;

S6. upgrade

and λ ^k+1, select new penalty coefficient

and μ ^k+1>=μ ^k; Parameter after upgrading is returned to step S4 to circulate.

In order to complete the reconstruction of PET image, the basic model of PET testing process is based on following equation:

y=Dx

Wherein, y is coincidence counting, here the data of measuring gained are generally sinogram, x is the values for spatial distribution of radioactive concentration in PET image, D is that (D is that m × n ties up matrix to system matrix, characterized transmitting photon and be detected the probability that device receives, it is subject to the decision of the factors such as panel detector structure, detection efficiency, decay, dead time).The process of rebuilding is collecting after y, calculates x by known sparse D.

Due to the sparse property of D, can not, directly by having inverted reconstruction, therefore consider transition problem model.

Further, consider that TV norm has good performance in image restoration, therefore set up the problem model based on TV norm, solve the problems referred to above by TV model, model is as follows:

$\begin{array}{cc}\underset{x}{\min }\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|\right|H_i\left(x\right)\left|\right|& s.t.\mathrm{Dx}=y\end{array}$

Wherein: H _i(x) be the discrete gradient on each pixel i of x, it is respectively x for two element values in bivector and this bivector _{i ,+1x}-x _iand x _{i ,+1y}-x _i; x _ifor the concentration value of i element in PET CONCENTRATION DISTRIBUTION vector x, x _{i ,+1x}for x in PET image corresponding to PET CONCENTRATION DISTRIBUTION vector x _ithe inner x that is expert at of institute _ithe concentration value of a rear column element, x _{i ,+1y}for x in PET image corresponding to PET CONCENTRATION DISTRIBUTION vector x _ix in column _ithe concentration value of a rear row element.

Due to the non-differentiability of the TV item of this model and non-linear, make the not recognition of this problem model again.Proceed the conversion of equivalence problem model.

Further, the enhancement mode Lagrange's equation that above formula is corresponding is:

$L_A(\mathrm{\&omega;},x)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\right|\left|{\mathrm{\&omega;}}_i\right||-v_i^T(H_i\left(x\right)-{\mathrm{\&omega;}}_i)+\frac{{\mathrm{\&beta;}}_i}{2}{\left|\right|H_i\left(x\right)-{\mathrm{\&omega;}}_i\left|\right|}_2^2)-{\mathrm{\&lambda;}}^T(\mathrm{Dx}-y)+\frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{Dx}-y\left|\right|}_2^2$

In conjunction with Fig. 2, in order to realize the correction slack variable in S5

with the image x rebuilding ^k+1, use alternating minimization method, comprise the following steps:

S501. initialization and setting

and x ⁰for initial point;

S502. calculate this renewal by following formula

${\mathrm{\&omega;}}_i^{k+1}=\max \left\{\right||H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\&beta;}}_i^k}||-\frac{1}{{\mathrm{\&beta;}}_i^k},0\}*\frac{(H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\&beta;}}_i^k})}{\left|\right|H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\&beta;}}_i^k}\left|\right|}$

S503. calculate the α of this renewal by following formula ^k, v ^kand λ ^k:

${\mathrm{\&alpha;}}^k=\frac{{(x^k-x^{k-1})}^T(x^k-x^{k-1})}{{(x^k-x^{k-1})}^T(Z^k\left(x^k\right)-Z^k\left(x^{k-1}\right))}$

$v_i^k=v_i^{k-1}-{\mathrm{\&beta;}}_i^{k-1}(H_i\left(x^{k-1}\right)-{\mathrm{\&omega;}}_i^{k-1})$

λ ^k=λ ^k-1-μ ^k-1(Dx ^k-1-y)

$Z^k\left(x^k\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[{\mathrm{\&beta;}}_i^k{H}_i^T(-H_i\left(x^k\right)-{\mathrm{\&omega;}}_i^{k+1})-H_i^T\left(v_i^k\right)]+{\mathrm{\&mu;}}^k{D}^T({\mathrm{Dx}}^k-y)-D^T{\mathrm{\&lambda;}}^k$

$Z^k\left(x^{k-1}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[{\mathrm{\&beta;}}_i^k{H}_i^T(-H_i\left(x^{k-1}\right)-{\mathrm{\&omega;}}_i^k)-H_i^T\left(v_i^k\right)]+{\mathrm{\&mu;}}^k{D}^T({\mathrm{Dx}}^{k-1}-y)-D^T{\mathrm{\&lambda;}}^k$

Wherein: for two element values in bivector and this bivector are respectively h _{i ,-1x}-h _iand h _{i ,-1y}-h _i; Wherein, h _ifor i element value of vector in (), h _{i ,-1x}for h in the PET image that in (), vector is corresponding _ithe inner h that is expert at of institute _iprevious column element value, h _{i ,-1y}for h in the PET image that in (), vector is corresponding _ih in column _iprevious row element value.

S504. calculate the x of this renewal by following formula _k+1:

x ^k+1=x ^k-α ^kZ ^k(x ^k)

If S505. the gap of twice renewal of judgement is enough little, no longer carry out next iteration, || x ^k+1-x ^k|| ₂enough little is the iteration stopping condition of alternating minimization, that is:

$\frac{{\left|\right|x^{k+1}-x^k\left|\right|}_2}{\left|\right|x^k\left|\right|}<\mathrm{\&rho;}$

Wherein: ρ is given convergence threshold, ρ=0.001 in present embodiment.

We adopt the validity of brain Voxel Phantom model experiment checking present embodiment below, and this model comprises some area with high mercury.This experiment running environment is: 8G internal memory, and 3.40GHz, 64 bit manipulation systems, CPU is intel Duo double-core.

PET method for reconstructing by present embodiment based on TV norm and traditional EM method reconstructed results compare, the two uses identical observed reading y and identical system matrix D to guarantee the comparability of result, design parameter arranges as follows: y is that n × n ties up the sinogram collecting, make m=n × n, D is that m × m ties up system matrix computed in advance, wherein β=2 ^-6, μ=2 ¹; Here n=64, i.e. m=4096.

For the checking of reconstructed image quality, adopt 64 projection angles of high-dimensional acquired original data, under each angle, beam is 64, i.e. m=4096, rebuilding image size is 64 × 64, i.e. dimension n=4096; Initial value is set the same.Image and traditional E M(greatest hope algorithm that Fig. 3 is true value image, rebuild based on present embodiment) the comparison schematic diagram of image rebuild of method, the image can visually see based on present embodiment reconstruction, compared with the result of EM, can recover more structure.In order further to study the degree of accuracy of rebuilding image, the image that the image that present embodiment is rebuild and traditional E M algorithm are rebuild has carried out the analysis of data, gets 32 row queue team data analysis, and data result as shown in Figure 4.

Table 1

	Present embodiment	EM
			Counting rate 1	19.67%	22.90%
Counting rate 2	23.34%	23.15%

For two kinds of different counting rate data, adopt respectively present embodiment and traditional EM method to compare, as shown in table 1; Application present embodiment reconstructed results with the overall relative error of true value on can reach essentially identical level with traditional EM method.Here counting rate 1 refers to and uses conventional normal counting rate data, is the 10^7 order of magnitude; Counting rate 2 refers to the data that use lower than conventional normal counting rate, is the 10^6 order of magnitude.As can be seen from Table 1, present embodiment, compared with traditional EM algorithm data reconstruction, in the situation that data volume declines, also can keep consistency substantially, and the feasibility to the low counting rate data reconstruction of technical solution of the present invention is described.

Claims (8)

1. the PET image rebuilding method based on TV norm, comprises the steps:

(1) utilize detector to survey the biological tissue that is injected with radiomaterial, collect the coincidence counting vector of current time, and this coincidence counting vector is proofreaied and correct;

(2), according to PET image-forming principle, the measurement equation of setting up PET is as follows:

y=Dx

Wherein: D is system matrix, y is coincidence counting after proofreading and correct and for m dimensional vector, and x is PET CONCENTRATION DISTRIBUTION vector and is n-dimensional vector, and m and n are the natural number that is greater than 1;

(3) introduce TV norm by the measurement equation to described, the TV problem model that obtains PET is as follows:

$\begin{array}{cc}\underset{x}{\min }\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|\right|H_i\left(x\right)\left|\right|& s.t.y=\mathrm{Dx}\end{array}$

Wherein: H _i(x) for two element values in bivector and this bivector are respectively x _{i ,+1x}-x _iand x _{i ,+1y}-x _i; x _ifor the concentration value of i element in PET CONCENTRATION DISTRIBUTION vector x, x _{i ,+1x}for x in PET image corresponding to PET CONCENTRATION DISTRIBUTION vector x _ithe inner x that is expert at of institute _ithe concentration value of a rear column element, x _{i ,+1y}for x in PET image corresponding to PET CONCENTRATION DISTRIBUTION vector x _ix in column _ithe concentration value of a rear row element, i is natural number and 1≤i≤n;

(4) described TV problem model is solved, the enhancement mode Lagrangian function that obtains TV problem model is as follows:

$L_A(\mathrm{\&omega;},x)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\right|\left|{\mathrm{\&omega;}}_i\right||-v_i^T(H_i\left(x\right)-{\mathrm{\&omega;}}_i)+\frac{{\mathrm{\&beta;}}_i}{2}{\left|\right|H_i\left(x\right)-{\mathrm{\&omega;}}_i\left|\right|}_2^2)-{\mathrm{\&lambda;}}^T(\mathrm{Dx}-y)+\frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{Dx}-y\left|\right|}_2^2$

Wherein: L _a(ω, x) is the enhancement mode Lagrangian function about ω and x, and ω is the gradient matrix corresponding with x, ω _ifor the two-dimensional gradient vector of corresponding i element in gradient matrix ω, v _iand β _ibe respectively the punishment vector sum penalty coefficient of corresponding i element, λ and μ are respectively the punishment vector sum penalty coefficient of Lagrangian function;

(5) described enhancement mode Lagrangian function is minimized and solved, estimate to obtain PET CONCENTRATION DISTRIBUTION vector x; And then carry out PET imaging according to the PET CONCENTRATION DISTRIBUTION vector x obtaining.

2. PET image rebuilding method according to claim 1, is characterized in that: in described step (4), utilize TVAL3 algorithm to solve TV problem model, obtain the enhancement mode Lagrangian function of TV problem model.

3. PET image rebuilding method according to claim 1, is characterized in that: in described step (5), by following iterative equation group, enhancement mode Lagrangian function is minimized and solved:

x ^k+1=x ^k-α ^kZ ^k(x ^k)

$\begin{array}{cc}{\mathrm{\&omega;}}_i^{k+1}=\max \left\{\right||H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\&beta;}}_i^k}||-\frac{1}{{\mathrm{\&beta;}}_i^k},0\}*\frac{(H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\&beta;}}_i^k})}{\left|\right|H_i\left(x^k\right)-\frac{v_i^k}{{\mathrm{\&beta;}}_i^k}\left|\right|}& {\mathrm{\&beta;}}_i^k\&GreaterEqual;{\mathrm{\&beta;}}_i^{k-1}\end{array}$

Wherein: x ^kand x ^k+1be respectively the PET CONCENTRATION DISTRIBUTION vector of the k time iteration and the k+1 time iteration,

be the gradient matrix ω of the k+1 time iteration ^k+1the two-dimensional gradient vector of middle corresponding i element,

with

be respectively the penalty coefficient of corresponding i element in the k time iteration and the k-1 time iteration,

be the punishment vector of corresponding i element in the k time iteration, α ^kbe the middle transition vector of the k time iteration, Z ^k(x ^k) be about x ^kpenalty term function; H _i() is that two element values in bivector and this bivector are respectively h _{i ,+1x}-h _iand h _{i ,+1y}-h _i; Wherein, h _ifor i element value of vector in (), h _{i ,+1x}for h in the PET image that in (), vector is corresponding _ithe inner h that is expert at of institute _ia rear column element value, h _{i ,+1y}for h in the PET image that in (), vector is corresponding _ih in column _irear a line element value.

4. PET image rebuilding method according to claim 3, is characterized in that: described middle transition vector α ^kcalculation expression as follows:

${\mathrm{\&alpha;}}^k=\frac{{(x^k-x^{k-1})}^T(x^k-x^{k-1})}{{(x^k-x^{k-1})}^T(Z^k\left(x^k\right)-Z^k\left(x^{k-1}\right))}$

Wherein: Z ^k(x ^k-1) be about x ^k-1penalty term function, x ^k-1it is the PET CONCENTRATION DISTRIBUTION vector of the k-1 time iteration.

5. PET image rebuilding method according to claim 4, is characterized in that: described penalty term function Z ^k(x ^k) and Z ^k(x ^k-1) calculation expression as follows:

$Z^k\left(x^k\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[{\mathrm{\&beta;}}_i^k{H}_i^T(-H_i\left(x^k\right)-{\mathrm{\&omega;}}_i^{k+1})-H_i^T\left(v_i^k\right)]+{\mathrm{\&mu;}}^k{D}^T({\mathrm{Dx}}^k-y)-D^T{\mathrm{\&lambda;}}^k$

$Z^k\left(x^{k-1}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[{\mathrm{\&beta;}}_i^k{H}_i^T(-H_i\left(x^{k-1}\right)-{\mathrm{\&omega;}}_i^k)-H_i^T\left(v_i^k\right)]+{\mathrm{\&mu;}}^k{D}^T({\mathrm{Dx}}^{k-1}-y)-D^T{\mathrm{\&lambda;}}^k$

Wherein:

be the gradient matrix ω of the k time iteration ^kthe two-dimensional gradient vector of middle corresponding i element, λ ^kand μ ^kbe respectively the punishment vector sum penalty coefficient of Lagrangian function in the k time iteration;

for two element values in bivector and this bivector are respectively h _{i ,-1x}-h _iand h _{i ,-1y}-h _i; Wherein, h _ifor i element value of vector in (), h _{i ,-1x}for h in the PET image that in (), vector is corresponding _ithe inner h that is expert at of institute _iprevious column element value, h _{i ,-1y}for h in the PET image that in (), vector is corresponding _ih in column _iprevious row element value.

6. PET image rebuilding method according to claim 3, is characterized in that: described punishment vector

calculation expression as follows:

$v_i^k=v_i^{k-1}-{\mathrm{\&beta;}}_i^{k-1}(H_i\left(x^{k-1}\right)-{\mathrm{\&omega;}}_i^{k-1})$

Wherein:

be the punishment vector of corresponding i element in the k-1 time iteration, x ^k-1be the PET CONCENTRATION DISTRIBUTION vector of the k-1 time iteration,

be the gradient matrix ω of the k-1 time iteration ^k-1the two-dimensional gradient vector of middle corresponding i element.

7. PET image rebuilding method according to claim 5, is characterized in that: described punishment vector λ ^kcalculation expression as follows:

λ ^k=λ ^k-1-μ ^k-1(Dx ^k-1-y) μ ^k≥μ ^k-1

Wherein: λ ^k-1and μ ^k-1be respectively the punishment vector sum penalty coefficient of Lagrangian function in the k-1 time iteration.

8. PET image rebuilding method according to claim 3, it is characterized in that: carry out iterative computation by described iterative equation group, the PET CONCENTRATION DISTRIBUTION vector x obtaining as estimation when the PET CONCENTRATION DISTRIBUTION vector reaching after maximum iteration time or iteration convergence; Iteration convergence condition is as follows:

$\frac{{\left|\right|x^{k+1}-x^k\left|\right|}_2}{\left|\right|x^k\left|\right|}<\mathrm{\&rho;}$

Wherein: ρ is given convergence threshold.

Images