CN107146218B - A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation - Google Patents
A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation Download PDFInfo
- Publication number
- CN107146218B CN107146218B CN201710233729.7A CN201710233729A CN107146218B CN 107146218 B CN107146218 B CN 107146218B CN 201710233729 A CN201710233729 A CN 201710233729A CN 107146218 B CN107146218 B CN 107146218B
- Authority
- CN
- China
- Prior art keywords
- matrix
- pet
- parameter
- rebuild
- tracer kinetics
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 44
- 239000000700 radioactive tracer Substances 0.000 title claims abstract description 38
- 238000003709 image segmentation Methods 0.000 title claims abstract description 11
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 15
- 230000011218 segmentation Effects 0.000 claims abstract description 15
- 239000011159 matrix material Substances 0.000 claims description 79
- 238000005259 measurement Methods 0.000 claims description 17
- 238000005457 optimization Methods 0.000 claims description 14
- 210000001519 tissue Anatomy 0.000 claims description 11
- 230000001360 synchronised effect Effects 0.000 claims description 6
- 230000017105 transposition Effects 0.000 claims description 6
- 239000012217 radiopharmaceutical Substances 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 3
- 210000002381 plasma Anatomy 0.000 claims description 2
- 238000004800 variational method Methods 0.000 claims 1
- 238000002474 experimental method Methods 0.000 abstract description 5
- 238000012636 positron electron tomography Methods 0.000 description 55
- 210000000115 thoracic cavity Anatomy 0.000 description 16
- 238000000342 Monte Carlo simulation Methods 0.000 description 13
- 230000006870 function Effects 0.000 description 11
- 230000002285 radioactive effect Effects 0.000 description 7
- 238000003384 imaging method Methods 0.000 description 5
- 238000012879 PET imaging Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 241000208340 Araliaceae Species 0.000 description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 235000008434 ginseng Nutrition 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000002600 positron emission tomography Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000000638 solvent extraction Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000003416 augmentation Effects 0.000 description 1
- 210000004027 cell Anatomy 0.000 description 1
- 210000001072 colon Anatomy 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 230000006806 disease prevention Effects 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000013399 early diagnosis Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 208000030159 metabolic disease Diseases 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10104—Positron emission tomography [PET]
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Nuclear Medicine (AREA)
Abstract
The dynamic PET images that the present invention provides a kind of based on image segmentation are rebuild and tracer kinetics method for parameter estimation, can obtain more accurately combining reconstructed results by combining statistics reconstruction model with the model for having physiological significance.When rebuilding and segmentation is coupled in a joint or solution frame simultaneously, for segmentation task, the information that noise model models Raw projection data can be obtained, and the result based on segmentation can also enhance the uniformity in each region, to realize the reconstructed results for more meeting truth.Compared with other individually rebuild the algorithm of dynamic PET images or estimated driving force parameter, the present invention can also obtain preferable reconstructed results.In conjunction with performance of the present invention in analogue data and truthful data experiment, compared with other individually rebuild the algorithm of dynamic PET images, estimated driving force parameter and image segmentation, the present invention can also obtain preferable reconstructed results.
Description
Technical field
The invention belongs to PET technical field of imaging, and in particular to a kind of dynamic PET images based on image segmentation rebuild and
Tracer kinetics method for parameter estimation.
Background technique
Positron emission tomography (Positron Emission Tomography, abbreviation PET) is nuclear medicine
One kind, its important role is just gradually shown in biomedical research and clinic diagnosis.PET imaging technique passes through to note
The increased radioactivity for entering tracer in organism is imaged, and the metabolic disorder of cell can be found from molecular level, is disease
The early diagnosis and prevention of disease provide effective foundation.Relative to radioactivity in bio-tissue can only be provided in certain time window
For the static PET imaging technique of intensity being evenly distributed, dynamic pet imaging, which is based on kinetic model, can obtain description biology
The quantitative function parameter of body vital movement has significant application value in scientific research and clinical application.
In the Problems of Reconstruction of dynamic PET, it is common practice to divide the image into, image reconstruction and Chemical kinetic parameter estimation
The problem independent as three separately solves.Requirement with dynamic pet imaging technology to temporal resolution is higher and higher, each
Limited photon counting constantly challenges dynamic pet imaging quality in frame measurement data.At this point, due to the photon of frame data
Count deficiency can not reflected well statistical property, each frame data are carried out using the algorithm for reconstructing based on statistical property single
The method solely rebuild often is unable to get very accurate reconstructed results.At this point, by combining tracer dynamics model to introduce
The prior information of radioactive intensity distribution from the time longitudinal axis, can be on the basis of improving PET reconstructed image quality simultaneously
Realize the estimation to kinetic parameter.In addition, image segmentation is also a conventional difficulties in PET imaging technique.In addition to right
It is the time-acttivity curve (time-activity to dynamic PET there are also a kind of way other than PET image is split
Curve, TAC) come the operation that is clustered, to distinguish functional area different in biological tissue.But either from PET's
On the plane of delineation still from the perspective of kinetic parameter to this problem, the partitioning algorithm accuracy based on reconstructed results is begun
The accuracy with back reconstructed results is relied on eventually, and partitioning algorithm can not be solved the problems, such as to measurement noise-sensitive.
In fact, the relationship between these three problems of image segmentation, image reconstruction and Chemical kinetic parameter estimation is very tight
Close.By the way that three is coupled into the same target equation, by selecting suitable noise model that can preferably match PET
Measurement data, segmentation result will reduce the susceptibility of noise at this time；And the introducing of segmentation result can make existing figure
Picture estimated result is more accurate, to generally obtain the reconstructed results for being more in line with truth.
Summary of the invention
In view of above-mentioned, the present invention provides a kind of dynamic PET images reconstruction and tracer kinetics ginseng based on image segmentation
Number estimation method, the segmentation and dynamic PET joint that can be achieved at the same time functional area are rebuild.
A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation, including walk as follows
It is rapid:
(1) biological tissue for being injected with radiopharmaceutical agent is detected using detector, dynamic acquisition obtains corresponding each
The coincidence counting vector at a moment, and establish coincidence counting matrix Y；
(2) dynamic PET image combined sequence is made to establish PET according to PET image-forming principle at PET concentration distribution matrix X
Measure equation；
(3) full variation (Total Variation, TV) constraint is introduced by measuring equation to PET, obtained based on TV's
PET image reconstruction model L (X)；
(4) using compartment model matching estimation tracer kinetics parameter, establish about X reconstruction model S synchronous with Ф's (X,
Ф)；
(5) region segmentation result obtained based on pretreatment, clusters tracer kinetics parameter matrix Ф, is gathered
Class parted pattern C (Ф)；
(6) it combines above three model L (X), S (X, Ф) to obtain the objective function LSC of synchronous reconstruction with C (Ф)
(X, Ф) is as follows:
Wherein: γ and ∈ is weight coefficient；
(7) PET concentration distribution matrix X is obtained after carrying out optimization to objective function LSC (X, Ф) and tracer is dynamic
Mechanics parameter matrix Ф.
The coincidence counting matrix Y is chronologically rearranged by each coincidence counting vector, the PET concentration distribution matrix X
It is chronologically rearranged by corresponding PET concentration distribution vector of each moment (i.e. a frame PET image).
The expression formula of the PET measurement equation is as follows:
Y=GX+R+S
Wherein: G is sytem matrix, and R and S are respectively the measurement noise matrix for reflecting chance event and scattering events.
The expression formula of the PET image reconstruction model L (X) is as follows:
Wherein: α is weight coefficient, and TV (X) is the full variation regular terms about X, gijIt is arranged for the i-th row jth in sytem matrix G
Element value (its indicate be the probability that is received by i-th of detector of photon being emitted at j-th of pixel), yimTo meet meter
I-th row m column element value, x in matrix number YjmFor jth row m column element value, r in PET concentration distribution matrix XimFor reflection with
I-th row m column element value, s in the measurement noise matrix R of machine eventimI-th in measurement noise matrix S to reflect scattering events
Row m column element value, i, j and m are natural number and 1≤i≤N, 1≤j≤K, 1≤m≤M, N are the dimension of coincidence counting vector
Degree, K are line number, that is, PET image pixel number of PET concentration distribution matrix X, and M is that the columns of PET concentration distribution matrix X is
Sampling time length.
The expression formula of the full variation regular terms TV (X) is as follows:
Wherein: D (xjm) it is about xjmTwo dimensional difference vector, the vector the first row element value be xjm-xJ, m+1, the second row
Element value is xjm-xJ+1, m, xJ, m+1For jth row m+1 column element value, x in PET concentration distribution matrix XJ+1, mFor PET concentration distribution
+ 1 row m column element value of jth in matrix X, | | | |2Indicate 2 norms.
The expression formula of the synchronous reconstruction model S (X, Ф) is as follows:
Wherein: μ is weight coefficient, and Ψ is dictionary matrix,TIndicate transposition, | | | |2Indicate 2 norms, | | | |1Indicate 1 norm.
The expression formula of the cluster segmentation MODEL C (Ф) is as follows:
Wherein: ChFor the parameter vector set for belonging to h class in tracer kinetics parameter matrix Ф, φhFor parameter vector collection
Close ChIn any parameter vector, nhFor parameter vector set ChIn parameter vector number, h is that natural number and 1≤h≤H, H are
The class number of cluster,TIndicate transposition.
Using ADMM, (Alternating Direction Method of Multipliers is handed in the step (7)
For direction Multiplier Algorithm) combine soft-threshold (Soft-Thresholding) iteration optimization algorithms to objective function LSC (X, Ф) into
Row optimization；Wherein, ADMM is iterated Optimization Solution for PET concentration distribution matrix X, and soft-threshold is directed to tracer power
Parameter matrix Ф iteration optimization is learned to solve.
The soft-threshold iteration optimization algorithms are based on following equation and carry out linear process to tracer kinetics parameter matrix
Ф is iterated update:
Nh=Ih-Eh/nh
Wherein: EhThe n all formed by 1 for oneh×nhTie up matrix, IhFor one and matrix EhThe identical unit square of dimension
Battle array, μ is weight coefficient,TIndicate that transposition, Ψ are dictionary matrix, ChFor the ginseng for belonging to h class in tracer kinetics parameter matrix Ф
Number vector set, φhFor parameter vector set ChIn any parameter vector, nhFor parameter vector set ChIn parameter vector
Number, xhFor in PET concentration distribution matrix X with parameter vector φhA corresponding TAC, | | | |2Indicate 2 norms, | | | |1It indicates
1 norm, h is natural number and 1≤h≤H, H are the class number of cluster.
The expression formula of the dictionary matrix Ψ is as follows:
Wherein: CI(t) and CI(τ) is respectively the concentration value of t moment and τ moment radiopharmaceutical agent in blood plasma,
With
Respectively about m group coincidence counting vector acquisition at the beginning of and the end time, θcCorrespond to c-th of chamber tissue index
When the columns that the coefficient of function, m and c are natural number and 1≤m≤M, 1≤c≤Z, M are PET concentration distribution matrix X samples
Between length, Z is natural number greater than 1；θ1~θNValue be in section [θmin, θmax] in chosen by exponential interval,
θminAnd θmaxRespectively the bound threshold value of coefficient, t and τ indicate the time.
The present invention can be obtained more accurately by combining statistics reconstruction model with the model for having physiological significance
Joint reconstructed results.It, can for segmentation task when rebuilding and segmentation is coupled in a joint or solution frame simultaneously
To obtain the information that noise model models Raw projection data, and the result based on segmentation can also enhance each region
Interior uniformity, to realize the reconstructed results for more meeting truth.Dynamic PET images are individually rebuild with other or are estimated
The algorithm of meter kinetic parameter is compared, and the present invention can also obtain preferable reconstructed results.In conjunction with the present invention in analogue data and
The calculation of dynamic PET images, estimated driving force parameter and image segmentation is individually rebuild in performance in truthful data experiment with other
Method is compared, and the present invention can also obtain preferable reconstructed results.
Detailed description of the invention
Fig. 1 is the template image of the thoracic cavity Monte Carlo simulation Zubal data.
Fig. 2 (a) is the true picture of the 4th frame of the thoracic cavity Zubal data.
Fig. 2 (b) is that data counts rate is using ML-EM method under 5*10^4 to the thoracic cavity Monte Carlo simulation Zubal data
The 4th frame image result rebuild.
Fig. 2 (c) is that data counts rate is using ML-EM method under 1*10^5 to the thoracic cavity Monte Carlo simulation Zubal data
The 4th frame image result rebuild.
Fig. 2 (d) is that data counts rate is using the method for the present invention under 5*10^4 to the thoracic cavity Monte Carlo simulation Zubal data
The 4th frame image result rebuild.
Fig. 2 (e) is that data counts rate is using the method for the present invention under 1*10^5 to the thoracic cavity Monte Carlo simulation Zubal data
The 4th frame image result rebuild.
Fig. 3 (a) is the true picture of the 6th frame of the thoracic cavity Zubal data.
Fig. 3 (b) is that data counts rate is using ML-EM method under 5*10^4 to the thoracic cavity Monte Carlo simulation Zubal data
The 6th frame image result rebuild.
Fig. 3 (c) is that data counts rate is using ML-EM method under 1*10^5 to the thoracic cavity Monte Carlo simulation Zubal data
The 6th frame image result rebuild.
Fig. 3 (d) is that data counts rate is using the method for the present invention under 5*10^4 to the thoracic cavity Monte Carlo simulation Zubal data
The 6th frame image result rebuild.
Fig. 3 (e) is that data counts rate is using the method for the present invention under 1*10^5 to the thoracic cavity Monte Carlo simulation Zubal data
The 6th frame image result rebuild.
Fig. 4 (a) is cluster result schematic diagram of the Monte Carlo simulation Zubal thoracic cavity data under using the method for the present invention.
Fig. 4 (b) is cluster result signal of the Monte Carlo simulation Zubal thoracic cavity data under using k mean value classification method
Figure.
Fig. 4 (c) is Monte Carlo simulation Zubal thoracic cavity data using based on poly- under kinetic model profile classification method
Class result schematic diagram.
Specific embodiment
In order to more specifically describe the present invention, with reference to the accompanying drawing and specific embodiment is to technical solution of the present invention
It is described in detail.
The present invention is based on the dynamic PET and kinetic parameter of region segmentation to combine method for reconstructing, includes the following steps:
(1) according to Model Establishment the measurement data matrix Y and sytem matrix G of dynamic PET scan.
The scanning process of dynamic PET is marked off into a certain number of time frames, detector in each time frame as required
The measurement data matrix Y that can go out a dynamic PET according to the sequential build of time of the coincidence counting vector collected；And
The probability that the photon being emitted at each pixel is received by each detector is counted, to obtain sytem matrix G.
(2) according to measurement data label tracer compound positron radionuclide type and time frame period distribution,
Dictionary matrix Ψ based on the theoretical settling time basic function of double compartment models.
If analogue data, then index parameters θ is arranged with the half-life period of position nucleic according to the tracer of simulation and detector is swept
It retouches time interval and establishes corresponding dictionary matrix；If truthful data, then according to tracer used in practical operation with position nucleic
Corresponding dictionary matrix is established with experimental setup.
The radioactive intensity that can be solved in tissue according to compartment model changes with time, instant m- radioactivity function
Be input tissue TAC and δ (t) function and a class index function convolution, be all tissue compartments TAC it
With.Therefore, a PET tracer compartment model can be indicated with one group of order-1 linear equation, and each equation indicates a difference
Chamber in tracer radioactive intensity, it may be assumed that
s.t.ψ0(t)=CI(t)
Wherein: CT(t) radioactive intensity in t moment tissue is indicated.N indicates the sum of different chambers, ψc(t) c is indicated
The function that the corresponding radioactive intensity of kind chamber changes over time, φcIt is ψc(t) corresponding coefficient.θcIt is the finger of corresponding chamber
Parameter in number function.CI(t) radioactive intensity of t moment input is indicated.
The dynamic PET images for each frame rebuild can be regarded as in the time interval of this frame and previous frame to chamber tissue
When m- activity curve integral, it may be assumed that
Wherein, xjmIt is radioactive intensity corresponding to the pixel j of m frame in PET image sequence,
With
It is m respectively
At the beginning of frame and the end time,It is corresponding TAC at the pixel j of mode tissue.Based on this relationship, we are by one
It is as follows that the corresponding equation group of group TAC is extended to a dictionary matrix:
Wherein:
(3) it initializes, weight coefficient α, γ and μ is set, the step-length σ in soft-threshold algorithm, the number of iterations k set 0, and setting is most
Big the number of iterations kMAX。
Weight coefficient α is [27.5,28.5] left and right, this is that the term of reference provided according to original existing algorithm is adjusted
It is whole to obtain.For weight coefficient γ generally in the range of [1,5] left and right, this is what the result based on experiment obtained.Weight system
Number μ generally can be selected to obtain suitably value (according to bibliography Positron in the range of [0.01,0.15]
Emission Tomography Compartmental Models&colon；A Basis Pursuit Strategy for
Experimental result in Kinetic Modeling), when noise increases, algorithm cannot be determined accurately just from dictionary
True time-acttivity curve needs suitably to increase at this time the value of μ, enhances sparse constraint.
(4) enter algorithm iteration renewal process, for kth time iteration, fixed kinetic parameter ΦkObtain PET image
Update result Xk+1。
4.1, which extract in target component the only part including PET image X, obtains following sub- optimization problem:
Wherein: TV (X) is full variation regular terms, and α is weight coefficient, and i=1 ..., N are the number of detector, m=1 ...,
M represents time frame, j=1 ..., and K represents the number of pixel on imaging plane；Therefore yimIt is the survey of i-th of detector of m frame
Measure data, xjmIt is the concentration value of j-th of pixel of m frame, and gijIt is then an element of sytem matrix, represents j-th of pixel
The probability that the photon being emitted at point is received by i-th of detector, rimAnd simRespectively indicate corresponding time frame and detector
Chance event and scattering events in measurement data.
4.2 calculation formula based on full variation regular terms define new variables ωjmAnd increasing is converted by the target equation in 4.1
Wide Lagrange's equation is as follows:
Wherein: wherein vector υjmIt is Lagrange multiplier, βjmIt is then penalty coefficient, xmIt is m frame PET image pixel group
At vector,
It is then xmDiscrete partially micro- operator at j-th of pixel.
4.3 fixed X, utilize two-dimensional contraction formula more new variables ωjm；
4.4 fixed ωjm, by LA(ωjm, X) and PET image X is updated according to following formula to after X derivation:
Xk+1=Xk-ρkdk(X)Xk
Wherein: ρkIt is the most fast decline step-length determined by reverse Non-monotone linear search, dkIt (X) is augmentation Lagrange
Derivative of the journey to X；
4.5 fixed X and ωjm, update Lagrange multiplier；
4.6 judge whether to meet iteration stopping condition: the variation < 10 of X-3Or reach maximum sub- the number of iterations, it is unsatisfactory for then
Return to 4.4；The iteration stopping if meeting, obtains the update result X of PET image Xk+1。
(5) for kth time iteration, fixed PET image Xk+1With kinetic parameter Φk, more using Di Li Cray cluster process
New cluster result
(6) for kth time iteration, fixed PET image Xk+1Obtain the estimated result Φ of kinetic parameterk+1:
6.1, which extract in target component the only part including kinetic parameter Φ, obtains following sub- optimization problem:
Wherein:Represent the subset for belonging to the pixel of h class in coefficient matrix Φ, nhIt is in the set
The number of pixel, ρhThe matrix being made of the average value of the corresponding pixel of h class, ∈ are the weight system for controlling submodel
Number.
6.2 by defining a new matrix Nh=Ih-Eh/nh, it is assumed that EhIt is the n all formed by 1h×nhDimension
Matrix, IhIt is one and EhThe identical unit matrix of dimension simultaneously defines new matrix Nh=Ih-Eh/nh.It introducesAbbreviation
Above-mentioned equation:
The kinetic parameter of the available update of soft-threshold iteration operator is utilized after 6.3 pairs of above formulas progress linearization process
Φk+1。
(7) judge whether to meet iteration stopping condition: the variation < 10 of X and Φ-3Or reach maximum number of iterations kMAX, no
Meet then return step (4), the iteration stopping if meeting obtains PET image X, tracer kinetics parameter Φ and cluster result.
We carry out experiment by the thoracic cavity the Zubal template data to Monte Carlo simulation to verify system of the present invention below
System is rebuild and the accuracy of segmentation result, and Fig. 1 is the template schematic diagram of the experiment thoracic cavity Zubal data used, by different areas
Domain is divided into three interested regions (region of interest, ROI).Test running environment are as follows: 8G memory, 3.40GHz,
64 bit manipulation systems, CPU are intel i7-3770；The PET scanner model Hamamatsu SHR-22000 simulated, if
Fixed radionuclide and drug be18F-FDG, setting sinogram are 128 projection angles, 128 beams under each angle
Collected data result, the size of sytem matrix G are 16384 × 16384.In this trial, to 5 × 104、1×105Two
Data for projection under the different counting rate of kind is tested.
To the weight of the result of PET image and traditional ML-EM (maximum likelihood-expectation maximum) in reconstruction framework of the present invention
It builds result to compare, the two is using identical measurement data matrix Y and sytem matrix G with the comparativity of control result.Fig. 2 (a)
It is distribution true value with Fig. 3 (a), it in counting rate is 5 × 10 that Fig. 2 (b)~Fig. 2 (e), which is in ML-EM respectively,4With 1 × 105The case where
Reconstruction framework lower and of the present invention is 5 × 10 in counting rate4With 1 × 105In the case where dynamic PET images sequence in the 4th frame
Reconstructed results；Fig. 3 (b)~Fig. 3 (e) be respectively ML-EM counting rate be 5 × 104With 1 × 105In the case where and the present invention
Reconstruction framework is 5 × 10 in counting rate4With 1 × 105In the case where dynamic PET images sequence in the 6th frame reconstructed results.With
Will become apparent from joint reclosing acquisition result in region noise be significantly less than ML-EM acquisition image, guarantee edge contrast
In the case where it is more smooth in functional area, table 1 is its further quantitative analysis result.
Table 1
Fig. 4 (a)~Fig. 4 (c) is the comparison result of cluster, the functional area distribution and have that reconstruction framework of the present invention obtains
The classification of k mean value and the obtained result of the profile classification method based on kinetic model be compared.It can be seen that joint is rebuild
Cluster result in noise it is significantly greater, this is mainly due to Di Li Cray process carry out classification pretreatment when be with pixel
Point is what unit was handled, therefore very sensitive to noise.
The above-mentioned description to embodiment is for that can understand and apply the invention convenient for those skilled in the art.
Person skilled in the art obviously easily can make various modifications to above-described embodiment, and described herein general
Principle is applied in other embodiments without having to go through creative labor.Therefore, the present invention is not limited to the above embodiments, ability
Field technique personnel announcement according to the present invention, the improvement made for the present invention and modification all should be in protection scope of the present invention
Within.
Claims (10)
1. a kind of dynamic PET images based on image segmentation are rebuild and tracer kinetics method for parameter estimation, include the following steps:
(1) biological tissue for being injected with radiopharmaceutical agent is detected using detector, when dynamic acquisition obtains corresponding to each
The coincidence counting vector at quarter, and establish coincidence counting matrix Y；
(2) dynamic PET image combined sequence is made to establish PET measurement according to PET image-forming principle at PET concentration distribution matrix X
Equation；
(3) full variational methods are introduced by measuring equation to PET, obtains the PET image reconstruction model L (X) based on TV；
(4) it using compartment model matching estimation tracer kinetics parameter, establishes about X reconstruction model S (X, Ф) synchronous with Ф's；
(5) region segmentation result obtained based on pretreatment, clusters tracer kinetics parameter matrix Ф, obtains cluster point
Cut MODEL C (Ф)；
(6) above three model L (X), S (X, Ф) and C (Ф) are combined to obtain the objective function LSC (X, Ф) of synchronous reconstruction
It is as follows:
Wherein: γ and ∈ is weight coefficient；
(7) PET concentration distribution matrix X and tracer kinetics are obtained after carrying out optimization to objective function LSC (X, Ф)
Parameter matrix Ф.
2. dynamic PET images according to claim 1 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The coincidence counting matrix Y is chronologically rearranged by each coincidence counting vector, and the PET concentration distribution matrix X is by each moment
Corresponding PET concentration distribution vector chronologically rearranges.
3. dynamic PET images according to claim 1 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The expression formula of the PET measurement equation is as follows:
Y=GX+R+S
Wherein: G is sytem matrix, and R and S are respectively the measurement noise matrix for reflecting chance event and scattering events.
4. dynamic PET images according to claim 3 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The expression formula of the PET image reconstruction model L (X) is as follows:
Wherein: α is weight coefficient, and TV (X) is the full variation regular terms about X, gijFor the i-th row jth column element in sytem matrix G
Value, yimFor the i-th row m column element value, x in coincidence counting matrix YjmFor jth row m column element in PET concentration distribution matrix X
Value, rimI-th row m column element value in measurement noise matrix R to reflect chance event, simFor the measurement for reflecting scattering events
I-th row m column element value in noise matrix S, i, j and m are natural number and 1≤i≤N, 1≤j≤K, 1≤m≤M, N are to meet
The dimension of count vector, K are line number, that is, PET image pixel number of PET concentration distribution matrix X, and M is PET concentration distribution
The columns of matrix X, that is, sampling time length.
5. dynamic PET images according to claim 4 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The expression formula of the full variation regular terms TV (X) is as follows:
Wherein: D (xjm) it is about xjmTwo dimensional difference vector, the vector the first row element value be xjm-xJ, m+1, the second row element
Value is xjm-xJ+1, m, xJ, m+1For jth row m+1 column element value, x in PET concentration distribution matrix XJ+1, mFor PET concentration distribution matrix
+ 1 row m column element value of jth in X, | | | |2Indicate 2 norms.
6. dynamic PET images according to claim 1 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The expression formula of the synchronous reconstruction model S (X, Ф) is as follows:
Wherein: μ is weight coefficient, and Ψ is dictionary matrix,TIndicate transposition, | | | |2Indicate 2 norms, | | | |1Indicate 1 norm.
7. dynamic PET images according to claim 1 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The expression formula of the cluster segmentation MODEL C (Ф) is as follows:
Wherein: ChFor the parameter vector set for belonging to h class in tracer kinetics parameter matrix Ф, φhFor parameter vector set Ch
In any parameter vector, nhFor parameter vector set ChIn parameter vector number, h be natural number and 1≤h≤H, H be cluster
Class number,TIndicate transposition.
8. dynamic PET images according to claim 1 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
Optimization is carried out to objective function LSC (X, Ф) using ADMM combination soft-threshold iteration optimization algorithms in the step (7)；
Wherein, ADMM is iterated Optimization Solution for PET concentration distribution matrix X, and soft-threshold is directed to tracer kinetics parameter matrix Ф
Iteration optimization solves.
9. dynamic PET images according to claim 8 are rebuild and tracer kinetics method for parameter estimation, it is characterised in that:
The soft-threshold iteration optimization algorithms are based on following equation and carry out linear process to change to tracer kinetics parameter matrix Ф
In generation, updates:
Nh=Ih-Eh/nh
Wherein: EhThe n all formed by 1 for oneh×nhTie up matrix, IhFor one and matrix EhThe identical unit matrix of dimension, μ
For weight coefficient,TIndicate that transposition, Ψ are dictionary matrix, ChFor belong in tracer kinetics parameter matrix Ф the parameter of h class to
Duration set, φhFor parameter vector set ChIn any parameter vector, nhFor parameter vector set ChIn parameter vector number,
xhFor in PET concentration distribution matrix X with parameter vector φhA corresponding TAC, | | | |2Indicate 2 norms, | | | |1Indicate 1 model
Number, h is natural number and 1≤h≤H, H are the class number of cluster.
10. dynamic PET images according to claim 6 or 9 are rebuild and tracer kinetics method for parameter estimation, feature exist
In: the expression formula of the dictionary matrix Ψ is as follows:
Wherein: CI(t) and CI(τ) is respectively the concentration value of t moment and τ moment radiopharmaceutical agent in blood plasma,
WithRespectively
For about m group coincidence counting vector acquisition at the beginning of and the end time, θcCorrespond to c-th of chamber tissue exponential function
Coefficient, columns, that is, sampling time that m and c are natural number and 1≤m≤M, 1≤c≤Z, M are PET concentration distribution matrix X is long
Degree, Z are the natural number greater than 1；θ1~θNValue be in section [θmin, θmax] in chosen by exponential interval, θminWith
θmaxRespectively the bound threshold value of coefficient, t and τ indicate the time.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710233729.7A CN107146218B (en) | 2017-04-11 | 2017-04-11 | A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710233729.7A CN107146218B (en) | 2017-04-11 | 2017-04-11 | A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107146218A CN107146218A (en) | 2017-09-08 |
| CN107146218B true CN107146218B (en) | 2019-10-15 |
Family
ID=59773553
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710233729.7A Active CN107146218B (en) | 2017-04-11 | 2017-04-11 | A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107146218B (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108627272B (en) * | 2018-03-22 | 2020-04-24 | 北京航空航天大学 | Two-dimensional temperature distribution reconstruction method based on four-angle laser absorption spectrum |
| CN108932741B (en) * | 2018-06-14 | 2022-08-19 | 上海联影医疗科技股份有限公司 | Dynamic PET parameter imaging method, device, system and computer readable storage medium |
| CN112365479B (en) * | 2020-11-13 | 2023-07-25 | 上海联影医疗科技股份有限公司 | PET parameter image processing method, device, computer equipment and storage medium |
| EP4331490A4 (en) * | 2021-06-11 | 2024-04-24 | Shanghai United Imaging Healthcare Co., Ltd. | Parameter imaging system and method |
| CN116671946A (en) * | 2023-04-14 | 2023-09-01 | 佛山读图科技有限公司 | Method for reconstructing dynamic image based on SPECT dynamic acquisition data |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009057019A1 (en) * | 2007-10-31 | 2009-05-07 | Koninklijke Philips Electronics N. V. | Tracer kinetic models for acoustic contrast imaging applications using photo-acoustics or thermo-acoustics |
| CN104657950A (en) * | 2015-02-16 | 2015-05-27 | 浙江大学 | Dynamic PET (positron emission tomography) image reconstruction method based on Poisson TV |
| CN105894550A (en) * | 2016-03-31 | 2016-08-24 | 浙江大学 | Method for synchronously reconstructing dynamic PET image and tracer kinetic parameter on the basis of TV and sparse constraint |
| CN106204674A (en) * | 2016-06-29 | 2016-12-07 | 浙江大学 | The dynamic PET images method for reconstructing retrained based on structure dictionary and kinetic parameter dictionary joint sparse |
| CN106251297A (en) * | 2016-07-19 | 2016-12-21 | 四川大学 | A kind of estimation based on multiple image fuzzy core the rebuilding blind super-resolution algorithm of improvement |
-
2017
- 2017-04-11 CN CN201710233729.7A patent/CN107146218B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009057019A1 (en) * | 2007-10-31 | 2009-05-07 | Koninklijke Philips Electronics N. V. | Tracer kinetic models for acoustic contrast imaging applications using photo-acoustics or thermo-acoustics |
| CN104657950A (en) * | 2015-02-16 | 2015-05-27 | 浙江大学 | Dynamic PET (positron emission tomography) image reconstruction method based on Poisson TV |
| CN105894550A (en) * | 2016-03-31 | 2016-08-24 | 浙江大学 | Method for synchronously reconstructing dynamic PET image and tracer kinetic parameter on the basis of TV and sparse constraint |
| CN106204674A (en) * | 2016-06-29 | 2016-12-07 | 浙江大学 | The dynamic PET images method for reconstructing retrained based on structure dictionary and kinetic parameter dictionary joint sparse |
| CN106251297A (en) * | 2016-07-19 | 2016-12-21 | 四川大学 | A kind of estimation based on multiple image fuzzy core the rebuilding blind super-resolution algorithm of improvement |
Non-Patent Citations (2)
| Title |
|---|
| Segmentation-Based Regularization of Dynamic SPECT Reconstruction;T.Humphries et al.;《IEEE Nuclear Science Symposium Conference Record》;20091231;全文 * |
| 基于MRF先验的PET图像重建与动力学参数估计;占杰;《万方数据库》;20100119;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107146218A (en) | 2017-09-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107146218B (en) | A kind of dynamic PET images reconstruction and tracer kinetics method for parameter estimation based on image segmentation | |
| CN105894550B (en) | A kind of dynamic PET images and tracer kinetics parameter synchronization method for reconstructing based on TV and sparse constraint | |
| CN107133997B (en) | A kind of dual tracer PET method for reconstructing based on deep neural network | |
| US11445992B2 (en) | Deep-learning based separation method of a mixture of dual-tracer single-acquisition PET signals with equal half-lives | |
| CN106887025B (en) | A method of the mixing tracer dynamic PET concentration distributed image based on stack self-encoding encoder is rebuild | |
| Iriarte et al. | System models for PET statistical iterative reconstruction: A review | |
| CN104657950B (en) | Dynamic PET (positron emission tomography) image reconstruction method based on Poisson TV | |
| CN108550172B (en) | PET image reconstruction method based on non-local characteristics and total variation joint constraint | |
| CN109993808B (en) | Dynamic double-tracing PET reconstruction method based on DSN | |
| CN108986916B (en) | Dynamic PET image tracer agent dynamics macro-parameter estimation method based on stacked self-encoder | |
| CN109636869A (en) | The dynamic PET images method for reconstructing constrained based on non local full variation and low-rank | |
| CN100412877C (en) | Computer simulation method for visualized information of substance metabolism functions inside human body | |
| CN107146263B (en) | A kind of dynamic PET images method for reconstructing based on the constraint of tensor dictionary | |
| Pedemonte et al. | A machine learning method for fast and accurate characterization of depth-of-interaction gamma cameras | |
| CN110197516A (en) | A kind of TOF-PET scatter correction method based on deep learning | |
| Cavalcanti et al. | Unmixing dynamic PET images with variable specific binding kinetics | |
| Hsiao et al. | Joint-MAP Bayesian tomographic reconstruction with a gamma-mixture prior | |
| Li et al. | Direct estimation of left ventricular ejection fraction via a cardiac cycle feature learning architecture | |
| WO2024109762A1 (en) | Pet parameter determination method and apparatus, and device and storage medium | |
| Jiang et al. | Smoothing dynamic positron emission tomography time courses using functional principal components | |
| Ding et al. | Dynamic SPECT reconstruction with temporal edge correlation | |
| Romaszko et al. | Massive dimensionality reduction for the left ventricular mesh | |
| Tong et al. | Tracer kinetics guided dynamic PET reconstruction | |
| Tong et al. | A robust state-space kinetics-guided framework for dynamic PET image reconstruction | |
| Gu et al. | Image-domain bootstrapping of pet time-course data for assessment of uncertainty in complex regional summaries of mapped kinetics |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |