CN108765509A - A kind of fast image reconstruction method for linear imaging system - Google Patents
A kind of fast image reconstruction method for linear imaging system Download PDFInfo
- Publication number
- CN108765509A CN108765509A CN201810496663.5A CN201810496663A CN108765509A CN 108765509 A CN108765509 A CN 108765509A CN 201810496663 A CN201810496663 A CN 201810496663A CN 108765509 A CN108765509 A CN 108765509A
- Authority
- CN
- China
- Prior art keywords
- prior information
- projection
- vector
- hyperplane
- adjustment
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 57
- 238000003384 imaging method Methods 0.000 title claims abstract description 32
- 239000013598 vector Substances 0.000 claims abstract description 93
- 238000005457 optimization Methods 0.000 claims abstract description 49
- 238000012804 iterative process Methods 0.000 claims abstract description 8
- 230000001133 acceleration Effects 0.000 claims description 19
- 238000005259 measurement Methods 0.000 claims description 7
- 238000013178 mathematical model Methods 0.000 claims description 5
- 239000011159 matrix material Substances 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 abstract description 23
- 238000004364 calculation method Methods 0.000 abstract description 6
- 238000013170 computed tomography imaging Methods 0.000 description 5
- 238000002591 computed tomography Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000001419 dependent effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000009659 non-destructive testing Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000033772 system development Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The present invention discloses a kind of fast image reconstruction method for linear imaging system, in each round cycle projection operation of traditional algebraic reconstruction technique, it increases and adjustment is accelerated to be operated with prior information optimization several times, adjustment is accelerated to operate the distance for being reduced iteration solution vector to optimal solution vector every time, significantly improve algorithm the convergence speed, and the increased additional storage space of algorithm and calculation amount very little, it is negligible, and iterative process is independent of previous iteration variable, it is easy to bind directly the fast image reconstruction that prior information realizes linear imaging system in the case of incomplete projection, in addition, the method of the present invention is suitably introduced into relaxation parameter to noise circumstance and has derived respective formula, expand the scope of application of method.
Description
Technical field
The invention belongs to the technical fields such as biomedical imaging, non-destructive testing, image reconstruction, are used for more particularly to one kind
The fast image reconstruction method of linear imaging system.
Background technology
Image information is that information content includes a branch maximum, that content is most abundant in information engineering subject field, modern
Information in objective world is transformed into figure by the technologies such as imaging system integrated optical, electronics, computer, machinery in various ways
As information, the extreme enrichment visual world of the mankind is obtained in various fields such as medical diagnosis, non-destructive testing, detection, investigations
It obtained and was widely applied.Currently, the characteristics of modern imaging system development of new techniques, is mainly reflected in digitlization, multifunction, more
In the integrated application of dimensionization and informationization etc., function and the performance of imaging are largely dependent upon used imaging
Algorithm.Therefore, imaging algorithm is the key that the hot spot of imaging system and imaging technique research.
It for many linear imaging systems, mathematically can often be modeled with a linear model, wherein comprising to big
The solution of scale system of linear equations, adoptable imaging algorithm have analytic method and iterative solution method.With computed tomography
(CT) for being imaged, CT imaging systems at this stage mostly use filtered back-projection (FBP), its advantage is that calculation amount is small, rebuild
Speed is fast, and good reconstructed image quality can be obtained when data for projection is complete.But in actually detected, due to many objective
Often there is the case where being difficult to detect complete projection data in reason.When data for projection is insufficient, the image of FBP algorithms reconstruction
Quality degradation, it is difficult to meet practical application request.
From FBP be the analytic method of representative it is different, algebraic reconstruction technique (ART) be imaged at the beginning the stage just by image reconstruction
Problem is converted into solution (extensive) system of linear equations problem, can be chosen according to different requirements different in conjunction with prior information
Object function is iterated solution, and high quality graphic can be still rebuild when under the conditions of sparse or incomplete data for projection.But algebraically
The disadvantage of reconstruction method is computationally intensive, rebuilds that speed is slower, the time cost of algorithm be hinder it to promote it is important because
Element.
In recent years, the progress of computer technology is that image reconstruction creates new hardware environment, the skill of arithmetic reconstruction method
Art advantage more highlights.In order to make algebraic reconstruction algorithm be widely used in actual imaging system, it is badly in need of research and development and is simple and efficient
Quick algorithm for reconstructing realize quick accurately image in incomplete data for projection in conjunction with image prior information.
As the immediate prior art, a kind of patent that the present inventor has applied " quick algebraically applied to CT imagings
Method for reconstructing " (application number 201510988996.6) discloses on the basis of traditional algebraic reconstruction technique, and arbitrary selected section is super flat
Face carries out further acceleration adjustment operation on these hyperplane to the projection solution vector obtained by traditional algebraic reconstruction technique,
Make the updated solution vector of adjustment former projection solution vector on the line of corresponding vector on the preceding an iteration hyperplane most
It is excellent.This method significantly improves algorithm the convergence speed, obtains advantageous effect, but directly applying to incomplete data for projection
Image reconstruction when there are still some shortcomings：Acceleration adjustment operation in algorithm is same with preceding an iteration dependent on current vector
Vectorial relationship on one hyperplane, if executing priori in conjunction with prior information on the basis of the current vector that each iteration obtains
Advance data quality operates, then increased due to prior information optimization operation to destroy this dependence to vectorial disturbance or variation
Relationship leads to the subsequent iterative calculation failure of algorithm.Therefore, patent (application number 201510988996.6) institute applied is public
The method held can be subject to certain restrictions when other prior information Optimizing operators are combined, and need to do special processing, this
Being unfavorable for the algorithm neatly combines prior information to realize the fast image reconstruction in the case of incomplete data for projection.
Invention content
To solve the above-mentioned problems, the present invention provides a kind of fast image reconstruction method for linear imaging system, should
Method directly carries out that adjustment is repeatedly accelerated to operate with prior information optimization inside an iteration, only introduces an auxiliary vector
With auxiliary scalar, hardly additionally increase algorithm amount of storage and calculation amount, in iterative process independent of previous iteration variable,
It is easy to bind directly the fast image reconstruction in the case of prior information realization incomplete data, this method introduces relaxation parameter simultaneously
The respective formula for having derived the update of auxiliary scalar and acceleration adjustment operation, is suitable for noise circumstance.The technical side that the present invention uses
Case is：
A kind of fast image reconstruction method for linear imaging system initially sets up the discretization number of linear imaging system
Model and prior information mathematical model are learned, converts the imaging problem of linear imaging system to extensive Solving Linear
It is extensive linear with the iterative solution of prior information Optimizing operator then in conjunction with adjustment operator is accelerated with prior information optimization problem
Equation group optimizes prior information, high quality graphic is rebuild under the conditions of sparse or incomplete data for projection；The big rule solved
Mould system of linear equations is expressed as Ax=b, wherein x is that the vector of image to be asked indicates, is the vectors of N × 1, by two-dimentional or three-dimensional
Image is rearranged to obtain as one-dimensional, and A is M × N sytem matrixes, with vectorial aiThe transposition of the i-th row vectors of representing matrix A, b M
× 1 measurement data vector, uses biIndicate i-th of element of vector b；Its solution procedure includes the following steps：
Step 1：Initialization：Set initial solution vector x0, relaxation parameter ρ, γ, and added inside loop iteration every time
The whole number T with prior information optimization operation of velocity modulation, and number of iterations k=0 is enabled, current vector (i.e. currently estimation solution vector) x=
x0, assist scalar d=0, auxiliary vector z=x0；
Step 2：Iterative process：Following steps are iteratively repeated until meeting the condition of convergence to kth time：
(2.1) hyperplane projection sequence is set, s is usedlIndicate the serial number of first of super projection plane, wherein l=1,2,
... M, M indicate the number of hyperplane, and enable t=1, and it is excellent with prior information that t indicates that iteration-internal carries out accelerating to adjust for the t time
Change operation；
(2.2) it is directed to each hyperplane i=sl, l=1,2 ... M, current vector x is executed successively it is following operate into
Row update：
(2.2.1) executes current vector x the projection operation to hyperplane i, updates current vector x and is marked with corresponding auxiliary
Measure d：
[x,d]←proj(x,d,ρ,ai,bi)
Wherein, [x', d']=proj (x, d, ρ, ai,bi), indicate that projection is operated with the d updates of auxiliary scalar, " ← " is to assign
It is worth symbol；
(2.2.2) ifAcceleration adjustment is then further executed to current vector x to operate with prior information optimization,
Update current vector x and corresponding auxiliary vector z and auxiliary scalar d：
[x,z,d]←lineAcc(x,z,d,γ)
t←t+1；
In above formula,Expression rounds up；[x', z', d']=lineAcc (x, z, d, γ) indicates to accelerate adjustment and priori
Advance data quality operates；
(2.3)k←k+1。
Further, in the step 1, the value range of relaxation parameter ρ, γ be (0,1], accelerate adjustment and priori
The value range that Advance data quality operates number T is 1≤T≤M.
Further, being thrown using fixed projection sequence or accidental projection sequence setting hyperplane in the step (2.1)
Shadow sequence.
Further, setting hyperplane projection sequence according to fixed projection sequence in the step (2.1), every time repeatedly
The serial number i=s of a super projection planes of generation l (l=1,2 ... M)lValue be fixed, then each iterative step (2.2.2)
In the optimization operation of acceleration adjustment and prior information be directed to identical hyperplane serial number；Conversely, suitable according to accidental projection
Sequence sets hyperplane projection sequence, the serial number i=s of first of super projection plane of each iterationlIt is random, then each iteration step
Suddenly the acceleration adjustment in (2.2.2) is directed to different hyperplane serial numbers from prior information optimization operation.
Further, projection operation [x', d']=proj (x, d, ρ, a in the step (2.2.1)i,bi) definition is such as
Under：
[x', d']=proj (x, d, ρ, ai,bi)
λ=bi-<ai,x>(1)
As ρ=1, projection operation [x', d']=proj (x, d, λ, ai,bi) output x ' be current vector x to hyperplane
The subpoint of i.After having executed projection operation, meet：||x*-x||2-||x*-x'||2=d-d', wherein x*It is the true of equation group
Real solution, i.e., the intersection of all hyperplane.
Further, the acceleration adjustment in the step (2.2.2) is operated with prior information optimization
[x', z', d']=lineAcc (x, z, d, γ) is defined as follows：
[x', z', d']=lineAcc (x, z, d, γ)
U=x-z (4)
xa=x+ γ β u (6)
X'=Sup (xa) (7)
Z'=x'(8)
D'=0 (9)
Accelerate adjustment with prior information optimization operation [x', z', d']=lineAcc (x, z, d, γ) formula (4), (5),
(6) it executes and accelerates adjustment operation, when | | x*-z||2-||x*-x||2(formula is executed in algorithm iteration process accelerates adjustment to=0-d
With satisfaction always before prior information optimization operation) and when γ=1, x that formula (6) obtainsaIt is vector x and corresponding auxiliary vector z
On line x is really solved away from equation group*Nearest point；As 0 < γ≤1, xaIt is truer solution x from equation group than the vector x before update*
Distance closer to.Wherein, x*It is the intersection of all hyperplane, i.e., has for i=1,2 ..., M：<ai,x*>=bi。[x',z',
D']=lineAcc (x, z, d, γ) formula (7) execute prior information optimization operation, x'=Sup (xa) return x ' compare xaWith figure
The prior information of picture is closer.The prior information of image can be many times retouched with certain mathematic(al) manipulation there are many representation
It states.For example, a kind of effective and widely used mathematical description mode is full variation (Total variance, TV) model of image
Number, the TV norms of general pattern are smaller.When describing prior information using TV norms, the priori in the method for the present invention
Advance data quality operates x'=Sup (xa) be the TV norms for reducing image, that is, the x ' returned compares xaWith smaller image TV norms.
The formula (8) of [x', z', d']=lineAcc (x, z, d, γ), (9) have updated auxiliary vector z and auxiliary scalar d respectively.It executes
Accelerate adjustment with after prior information optimization operation [x', z', d']=lineAcc (x, z, d, γ), all things considered is updated
Vector x ' than the vector x before update has smaller measurement error, with the prior information of image more close to (such as TV norms are more
It is small).
The present invention theory analysis be：
The basic thought of algebraic reconstruction technique is：The hyperplane that each equation is regarded as to N-dimensional space, from initial solution vector
x0It sets out, in each iteration, according to certain projection sequence, carries out a wheel cycle projection operation to each hyperplane successively, gradually
Approach the true solution of equation group.By taking sequential projection (the hyperplane serial number 1,2 ..., the M that project successively) as an example, changing every time
Generation, by currently estimating solution vector (referred to as current vector) x, by hyperplane projection sequence i=1,2 ..., M are executed successively
Estimation solution vector x is constantly updated in following operation：
λ←bi-<ai,x> (10)
In above formula<·>For interior product code, " ← " is assignment.Formula (11) describes projection operation, on the right of formula
Vector x is the vector executed before projection operation, and the vector x on the formula left side is to execute newer vector after projection operation.ρ is ART
The relaxation parameter of algorithm, be applied to noise circumstance, value range generally (0,1] between.As ρ=1, formula (10), (11) will
Current vector x projects to the hyperplane of serial number i (on abbreviation hyperplane i).
In the methods of the invention, if the acceleration adjustment of step (2.2.2) and prior information optimization operation are removed, or will be into
Row accelerates adjustment and the number T of prior information optimization operation to be set as 0, then it is traditional algebraic reconstruction technique that the method for the present invention, which is degenerated,.
In addition after the optimization operation of the acceleration adjustment of step (2.2.2) and prior information, algorithm speed is improved, while sparse or not
High quality graphic is rebuild under the conditions of complete data for projection.
Illustrate first in step (2.2.1), the projection operation to hyperplane i is executed to current vector x
[x', d']=proj (x, d, ρ, ai,bi) after, have：
||x*-x||2-||x*-x'||2=d-d'(12)
It proves as follows：
As the relaxation parameter ρ=1 in formula (2), the vector x that is obtained by formula (2) ' (be expressed as xp) it is that current vector x exists
The subpoint of hyperplane i：
That is, xpMeet：
<|x-xp,x*-xp>=0 (14)
It is corresponding to have：
||x*-x||2=| | x*-xp||2+||x-xp||2 (15)
||x*-x'||2=| | x*-xp||2+||x'-xp||2 (16)
It can be obtained by formula (2) and formula (13)：
Convolution (15), (16), (17), (18) can obtain：
It can be obtained by formula (19) and formula (3)：
||x*-x||2-||x*-x'||2=d-d'
Proof finishes.
During algorithm iteration, vector x is corresponding with auxiliary scalar d, and auxiliary vector z is corresponding with auxiliary scalar 0.In conjunction with
Above-mentioned property, it is known that executed in step (2.2.2) and meet following conditions before accelerating adjustment and prior information optimization operation：
||x*-z||2-||x*-x||2=0-d (20)
Illustrate below, executed if γ=1, in step (2.2.2) accelerate adjustment and prior information optimization operation [x', z',
D']=lineAcc (x, z, d, γ) formula (6) after, the x of returnaVector x on the line of corresponding auxiliary vector z away from equation
The true solution x of group*Nearest point, that is, meet
<U, x*-xa>=0 (21)
U in above formula, xaRespectively by formula (4), (6) define (enabling γ=1), i.e., as γ=1, if the formula (6) by updating x
The x of returnaIt is that vector x really solves x on the line of corresponding auxiliary vector z away from equation group*Nearest point, then the β in formula (6) by
Formula (5) defines.
It proves as follows：
Convolution (21), (4), (6) can obtain：
||x*-x||2=| | x*-xa||2+β2||u||2 (22)
||x*-z||2=| | x*-xa||2+(1+β)2||u||2 (23)
By formula (22), (23), it is known that：
||x*-z||2-||x*-x||2=(2 β+1) | | u | |2 (24)
Convolution (20), (24), it is known that：
The as calculation expression of formula (5).
Proof finishes.
In view of noise reason, relaxation parameter ρ, γ are introduced in formula (2) and formula (6).As 0 < γ≤1, adjusted by accelerating
The whole vector x returned with prior information optimization operation [x', z', d']=lineAcc (x, z, d, γ) formula (6)aBefore update
Vector x really solves x from equation group*Distance closer to.
Using the advantageous effect of the technical program：
1. with the prior art (applying for a patent the method disclosed in (application number 201510988996.6) before the present inventor)
It compares, the method for the present invention directly carries out that adjustment is repeatedly accelerated to operate with prior information optimization in one cycle iteration-internal, due to
Iterative process is independent of previous iteration variable, it is easier to bind directly prior information realize it is fast in the case of incomplete projection
Fast image reconstruction.
2. with the prior art (applying for a patent the method disclosed in (application number 201510988996.6) before the present inventor)
It compares, the method for the present invention is suitably introduced into relaxation parameter to noise circumstance and has derived the update of auxiliary scalar and accelerated adjustment operation
Respective formula, expand the scope of application of method.
3. improving image reconstruction algorithm convergence rate：The method of the present invention merely adds one compared with traditional algebraic reconstruction technique
Auxiliary vector z, the increased additional storage space of institute can be ignored completely.In practical application, each iteration-internal accelerate adjustment with
The number T of prior information optimization operation can set to obtain very little (such as less than 10), and increased calculation amount additional in this way can equally neglect
Slightly, and simulation result shows that convergence speed of the algorithm can but significantly improve.
4. being capable of effectively extended parallel port system application：The method of the present invention can be directed to different linear imaging system applications, flexibly
Using corresponding prior information, the fast image reconstruction under the conditions of incomplete data is realized.
Description of the drawings
Invention is further described in detail below in conjunction with the accompanying drawings.
Fig. 1 is a kind of fast image reconstruction method flow schematic diagram for linear imaging system of the present invention；
Fig. 2 is that hyperplane i=s is directed in the method for the present inventionl, the flow diagram of update operation is executed to current vector x；
Fig. 3 is the comparison figure of one embodiment reconstruction image and original image of the method for the present invention；
Fig. 4 is the embodiment of the method for the present invention and combines the algorithm square vector error convergence of ART and TV norm optimizations bent
Line chart.
Specific implementation mode
To make the objectives, technical solutions, and advantages of the present invention clearer, right in the following with reference to the drawings and specific embodiments
The present invention is further elaborated.
In the present embodiment, referring to shown in Fig. 1, Fig. 2, a kind of fast image reconstruction method for CT imaging systems is first
The discrete mathematical model and prior information mathematical model for first establishing CT imaging systems, convert on a large scale CT imaging problems to
Solving Linear and prior information optimization problem, then in conjunction with acceleration adjustment operator and prior information Optimizing operator iteration
Extensive system of linear equations is solved, optimizes prior information, high-quality is rebuild under the conditions of sparse or incomplete data for projection
Picture, the extensive system of linear equations solved are expressed as Ax=b, wherein x be image to be asked vector indicate, be N × 1 to
Amount, by rearranging to obtain by one-dimensional to two dimensional image, A is M × N sytem matrixes, with vectorial aiThe i-th row vectors of representing matrix A
Transposition, b is M × 1 measurement data vector, uses biIndicate i-th of element of vector b；Its solution procedure includes the following steps：
Step 1：Initialization：Set initial solution vector x0, relaxation parameter ρ, γ, and added inside loop iteration every time
The whole number T with prior information optimization operation of velocity modulation, and enable number of iterations k=0, current vector x=x0, auxiliary scalar d=0 is auxiliary
Help vectorial z=x0；
In the step 1, the value range of relaxation parameter ρ, γ be (0,1], accelerate adjustment to be operated with prior information optimization
The value range of number T is 1≤T≤M.
Step 2：Iterative process：Following steps are iteratively repeated until meeting the condition of convergence to kth time：
(2.1) hyperplane projection sequence is set, s is usedlIndicate the serial number of first of super projection plane, wherein l=1,2,
... M, M indicate the number of hyperplane, and enable t=1, and t carries out accelerating adjustment and priori letter for the t time for indicating iteration-internal
Breath optimization operation；
Using fixed projection sequence or accidental projection sequence setting hyperplane projection sequence in the step (2.1).
(2.2) it is directed to each hyperplane i=sl, l=1,2 ... M, current vector x is executed successively it is following operate into
Row update：
(2.2.1) executes current vector x the projection operation to hyperplane i, updates current vector x and is marked with corresponding auxiliary
Measure d：
[x,d]←proj(x,d,ρ,ai,bi)
Wherein, [x', d']=proj (x, d, ρ, ai,bi), indicate that projection is operated with the d updates of auxiliary scalar, " ← " is to assign
It is worth symbol；
Projection operation [x', d']=proj (x, d, ρ, a in the step (2.2.1)i,bi) be defined as follows：
[x', d']=proj (x, d, ρ, ai,bi)
λ=bi-<ai,x> (1)
As ρ=1, projection operation [x', d']=proj (x, d, λ, ai,bi) output x ' be current vector x to hyperplane
The subpoint of i.After having executed projection operation, meet：||x*-x||2-||x*-x'||2=d-d', wherein x*It is the true of equation group
Real solution.
(2.2.2) ifAcceleration adjustment is then further executed to current vector x to grasp with prior information optimization
Make, updates current vector x and corresponding auxiliary vector z and auxiliary scalar d：
[x,z,d]←lineAcc(x,z,d,γ)
t←t+1；
In above formula,Expression rounds up；[x', z', d']=lineAcc (x, z, d, γ) indicates to accelerate adjustment and priori
Advance data quality operates；
Acceleration adjustment in the step (2.2.2) and prior information optimization operation [x', z', d']=lineAcc (x, z,
D, γ) it is defined as follows：
[x', z', d']=lineAcc (x, z, d, γ)
U=x-z (4)
xa=x+ γ β u (6)
X'=Sup (xa) (7)
Z'=x'(8)
D'=0 (9)
Accelerate adjustment with prior information optimization operation [x', z', d']=lineAcc (x, z, d, γ) formula (4), (5),
(6) it executes and accelerates adjustment operation, when | | x*-z||2-||x*-x||2(formula is executed in algorithm iteration process accelerates adjustment to=0-d
With satisfaction always before prior information optimization operation) and when γ=1, x that formula (6) obtainsaIt is vector x and corresponding auxiliary vector z
On line x is really solved away from equation group*Nearest point；As 0 < γ≤1, xaIt is truer solution x from equation group than the vector x before update*
Distance closer to.The formula (7) of [x', z', d']=lineAcc (x, z, d, γ) executes prior information optimization operation, x'=Sup
(xa) return x ' generally compare xaIt is closer with the prior information of image.There are many representations for the prior information of image, very much
When can be described with certain mathematic(al) manipulation, such as full variation (Total variance, TV) norm.It is retouched when using TV norms
When stating prior information, the prior information optimization operation x'=Sup (x in the method for the present inventiona) be the TV norms for reducing image, i.e.,
The x ' of return compares xaWith smaller image TV norms.The formula (8) of [x', z', d']=lineAcc (x, z, d, γ), (9) point
Auxiliary vector z and auxiliary scalar d are not had updated.Execute accelerate adjustment with prior information optimization operation [x', z', d']=
After lineAcc (x, z, d, γ), all things considered, updated vector x ' missed with smaller measurement than the vector x before update
Difference, the prior information with image is more close to (such as TV norms smaller).
The step (2.1) sets hyperplane projection sequence using fixed projection sequence, each iteration l (l=1,
2 ... M) a super projection plane serial number i=slValue be fixed, the then acceleration in each iteration below step (2.2.2)
Adjustment is directed to identical hyperplane serial number with prior information optimization operation；Conversely, being set according to accidental projection sequence super
Plane projection sequence, the serial number i=s of first of super projection plane of each iterationlIt is random, then each iteration below step
Acceleration adjustment in (2.2.2) is directed to different hyperplane serial numbers from prior information optimization operation.
(2.3) k ← k+1, i.e. number of iterations add 1, are ready for next iteration.
With reference to 1 specific example, the present invention is described further.
The present embodiment corresponds to the case where incomplete projection data, and CT images to be reconstructed are one 128 × 128 X-Y schemes
Picture, i.e. original image in Fig. 3, it is corresponding be expressed as 16384 × 1 vector x*, A is the system square that dimension is 7200 × 16384
Battle array, i.e. M=7200, N=16384, b=Ax*The measurement data vector that+e is 16384 × 1 (contains a small amount of complete zero in matrix A
Row, when processing will reject these data).E is zero mean Gaussian white noise vector, and amplitude intensity is the 1% of measurement data b.
In incomplete projection, in order to rebuild high quality graphic, it is also necessary to introduce some prioris and carry out canonical
Change.Full variation (TV) Regularization Technique is used in this example, CT image reconstructions problem is converted into following object function after regularization
Minimization problem：
In the present embodiment, the parameter of η=0.004, initial phase the method for the present invention is set as：ρ=1, γ=1, just
Beginning vector x0It is set as full null vector, the acceleration adjustment carried out in each loop iteration and prior information optimization operation number T=4,
Using accidental projection sequence.
In each iterative process, one group 1 is generated first, 2 ... the random alignment combination of M is used as projection sequence, then presses
The projection sequence randomly generated carries out projection to each hyperplane and is operated with prior information optimization with adjustment is accelerated.
(i.e. for the 1800th, 3600,5400,7200) other than super projection plane, only
Projection operation (2.2.1) is executed, acceleration adjustment and prior information optimization operation (2.2.2) are not executed；For the 1800th, 3600,
5400,7200 super projection plane also executes after executing projection operation (2.2.1) and adjustment is accelerated to be operated with prior information optimization
(2.2.2).Due to using accidental projection sequence, each iteration the 1800th, 3600,5400, the sequences of 7200 super projection planes
It is number general different.For example, in the serial number 353 of the 1800th super projection plane of first time iteration, in second of iteration the 1800th
The serial number 3478 of a super projection plane, this expression carry out acceleration adjustment in first time iteration to the hyperplane of serial number 353
Optimize with prior information and operate, and the hyperplane of serial number 3478 is carried out in second of iteration adjustment is accelerated with priori to believe
Breath optimization operation.
Due to the method for the present invention in iterative process independent of previous iteration variable, be easy to bind directly prior information
Realize the fast image reconstruction in the case of incomplete projection.The comparison of image and original image that the method for the present invention is rebuild
Figure as shown in figure 3, square vector error convergence curve figure of the method for the present invention and the algorithm for combining ART and TV norm optimizations such as
Shown in Fig. 4, as can be seen from the figure the method for the present invention significantly improves convergence rate.
The basic principles, main features and advantages of the present invention have been shown and described above.The technology of the industry
Personnel are it should be appreciated that the present invention is not limited to the above embodiments, and the above embodiments and description only describe this
The principle of invention, without departing from the spirit and scope of the present invention, various changes and improvements may be made to the invention, these
Changes and improvements all fall within the protetion scope of the claimed invention.The claimed scope of the invention is by appended claims
And its equivalent thereof.
Claims (6)
1. a kind of fast image reconstruction method for linear imaging system, which is characterized in that initially set up linear imaging system
Discrete mathematical model and prior information mathematical model, convert the imaging problem of linear imaging system to extensive linear side
Journey group solves and prior information optimization problem, is advised greatly then in conjunction with acceleration adjustment operator and the iterative solution of prior information Optimizing operator
Mould system of linear equations optimizes prior information, high quality graphic is rebuild under the conditions of sparse or incomplete data for projection；It is solved
Extensive system of linear equations is expressed as Ax=b, wherein x is that the vector of image to be asked indicates, is the vectors of N × 1, by two dimension or
3-D view is rearranged to obtain as one-dimensional, and A is M × N sytem matrixes, with vectorial aiThe transposition of the i-th row vectors of representing matrix A, b
For the measurement data vectors of M × 1, b is usediIndicate i-th of element of vector b；Its solution procedure includes the following steps：
Step 1：Initialization：Set initial solution vector x0, relaxation parameter ρ, γ, and every time iteration-internal carry out accelerate adjustment with first
The number T of Advance data quality operation is tested, and enables number of iterations k=0, current vector x=x0, assist scalar d=0, auxiliary vector z=x0；
Step 2：Iterative process：Following steps are iteratively repeated until meeting the condition of convergence to kth time：
(2.1) hyperplane projection sequence is set, s is usedlIndicate the serial number of first of super projection plane, wherein l=1,2 ... M,
M indicates the number of hyperplane, and enables t=1, and t indicates that iteration-internal carries out accelerating adjustment and prior information optimization for the t time
Operation；
(2.2) it is directed to each hyperplane i=sl, l=1,2 ... M executes current vector x following operation and carries out more successively
Newly：
(2.2.1) executes current vector x the projection operation to hyperplane i, updates current vector x and corresponding auxiliary scalar d：
[x,d]←proj(x,d,ρ,ai,bi)
Wherein, [x', d']=proj (x, d, ρ, ai,bi), indicate that projection is operated with the d updates of auxiliary scalar, " ← " is that assignment accords with
Number；
(2.2.2) ifAcceleration adjustment is then further executed to current vector x to operate with prior information optimization,
Update current vector x and corresponding auxiliary vector z and auxiliary scalar d：
[x,z,d]←lineAcc(x,z,d,γ)
t←t+1；
In above formula,Expression rounds up；[x', z', d']=lineAcc (x, z, d, γ) indicates to accelerate adjustment and prior information
Optimization operation；
(2.3)k←k+1。
2. a kind of fast image reconstruction method for linear imaging system according to claim 1, which is characterized in that institute
State in step 1, the value range of relaxation parameter ρ, γ be (0,1], accelerate adjustment to operate taking for number T with prior information optimization
It is worth ranging from 1≤T≤M.
3. a kind of fast image reconstruction method for linear imaging system according to claim 1, which is characterized in that institute
It states in step (2.1) using fixed projection sequence or accidental projection sequence setting hyperplane projection sequence.
4. a kind of fast image reconstruction method for linear imaging system according to claim 3, which is characterized in that institute
It states in step (2.1) and sets hyperplane projection sequence, the sequence of first of super projection plane of each iteration according to fixed projection sequence
Number i=slValue be fixed, then the acceleration adjustment in each iterative step (2.2.2) is directed to prior information optimization operation
It is identical hyperplane serial number；Conversely, setting hyperplane projection sequence, first of the throwing of each iteration according to accidental projection sequence
The serial number i=s of shadow hyperplanelIt is random, then the acceleration adjustment in each iterative step (2.2.2) is grasped with prior information optimization
It is directed to different hyperplane serial numbers.
5. a kind of fast image reconstruction method for linear imaging system according to claim 1, which is characterized in that institute
State the projection in step (2.2.1) and auxiliary scalar d update operation [x', d']=proj (x, d, ρ, ai,bi) be defined as follows：
[x', d']=proj (x, d, ρ, ai,bi)
λ=bi-<ai,x>
6. a kind of fast image reconstruction method for linear imaging system according to claim 1, which is characterized in that institute
The acceleration adjustment stated in step (2.2.2) is defined with prior information optimization operation [x', z', d']=lineAcc (x, z, d, γ)
It is as follows：
[x', z', d']=lineAcc (x, z, d, γ)
U=x-z
xa=x+ γ β u
X'=Sup (xa)
Z'=x'
D'=0.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810496663.5A CN108765509B (en) | 2018-05-22 | 2018-05-22 | Rapid image reconstruction method for linear imaging system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810496663.5A CN108765509B (en) | 2018-05-22 | 2018-05-22 | Rapid image reconstruction method for linear imaging system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108765509A true CN108765509A (en) | 2018-11-06 |
| CN108765509B CN108765509B (en) | 2020-08-14 |
Family
ID=64007545
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810496663.5A Expired - Fee Related CN108765509B (en) | 2018-05-22 | 2018-05-22 | Rapid image reconstruction method for linear imaging system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108765509B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115619890A (en) * | 2022-12-05 | 2023-01-17 | 山东省计算中心（国家超级计算济南中心） | Tomography method and system for solving linear equation set based on parallel random iteration |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090185636A1 (en) * | 2008-01-23 | 2009-07-23 | Sparsense, Inc. | Parallel and adaptive signal processing |
| US20140219529A1 (en) * | 2013-02-01 | 2014-08-07 | Toshiba Medical Systems Corporation | Alterative noise map estimation methods for ct images |
| CN105590332A (en) * | 2015-12-24 | 2016-05-18 | 电子科技大学 | Rapid algebraic reconstruction technique applied to computed tomography imaging |
| CN105608719A (en) * | 2015-12-28 | 2016-05-25 | 电子科技大学 | Rapid CT image reconstruction method based on two-stage projection adjustment |
| CN106097291A (en) * | 2016-06-08 | 2016-11-09 | 天津大学 | A kind of ART flame based on gradient total variation section restructing algorithm |
| CN106709879A (en) * | 2016-12-08 | 2017-05-24 | 中国人民解放军国防科学技术大学 | Spatial variation point diffusion function smoothing method based on simple lens calculating imaging |
-
2018
- 2018-05-22 CN CN201810496663.5A patent/CN108765509B/en not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090185636A1 (en) * | 2008-01-23 | 2009-07-23 | Sparsense, Inc. | Parallel and adaptive signal processing |
| US20140219529A1 (en) * | 2013-02-01 | 2014-08-07 | Toshiba Medical Systems Corporation | Alterative noise map estimation methods for ct images |
| CN105590332A (en) * | 2015-12-24 | 2016-05-18 | 电子科技大学 | Rapid algebraic reconstruction technique applied to computed tomography imaging |
| CN105608719A (en) * | 2015-12-28 | 2016-05-25 | 电子科技大学 | Rapid CT image reconstruction method based on two-stage projection adjustment |
| CN106097291A (en) * | 2016-06-08 | 2016-11-09 | 天津大学 | A kind of ART flame based on gradient total variation section restructing algorithm |
| CN106709879A (en) * | 2016-12-08 | 2017-05-24 | 中国人民解放军国防科学技术大学 | Spatial variation point diffusion function smoothing method based on simple lens calculating imaging |
Non-Patent Citations (2)
| Title |
|---|
| CHUAN LIN 等: "Algebraic Reconstruction Technique with adaptive relaxation parameter based on hyperplane distance and data noise level", 《2016 IEEE MTT-S INTERNATIONAL CONFERENCE ON NUMERICAL ELECTROMAGNETIC AND MULTIPHYSICS MODELING AND OPTIMIZATION (NEMO)》 * |
| 阙介民 等: "基于不完备投影数据重建的四种迭代算法比较研究", 《CT理论与应用研究》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115619890A (en) * | 2022-12-05 | 2023-01-17 | 山东省计算中心（国家超级计算济南中心） | Tomography method and system for solving linear equation set based on parallel random iteration |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108765509B (en) | 2020-08-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sun et al. | Block coordinate regularization by denoising | |
| Hawe et al. | Dense disparity maps from sparse disparity measurements | |
| Chen et al. | Adaptively regularized constrained total least-squares image restoration | |
| CN108765514B (en) | Acceleration method and device for CT image reconstruction | |
| CN110060315B (en) | Image motion artifact eliminating method and system based on artificial intelligence | |
| CN105488767B (en) | A kind of compressed sensing image fast reconstructing method based on Least-squares minimization | |
| CN111870245A (en) | Cross-contrast-guided ultra-fast nuclear magnetic resonance imaging deep learning method | |
| CN108765540B (en) | Relighting method based on image and ensemble learning | |
| CN112634149A (en) | Point cloud denoising method based on graph convolution network | |
| CN114913262B (en) | Nuclear magnetic resonance imaging method and system with combined optimization of sampling mode and reconstruction algorithm | |
| Jiang et al. | Cnn driven sparse multi-level b-spline image registration | |
| Shen et al. | 3D shape reconstruction from images in the frequency domain | |
| CN114283235B (en) | Three-dimensional magnetic layer reconstruction method and system based on limited angle projection data | |
| Dolmatova et al. | Accelerated FBP for computed tomography image reconstruction | |
| Malhotra et al. | Tomographic reconstruction from projections with unknown view angles exploiting moment-based relationships | |
| CN108765509A (en) | A kind of fast image reconstruction method for linear imaging system | |
| CN105608719B (en) | A kind of rapid CT image rebuilding method based on two benches projection adjustment | |
| CN110060314B (en) | CT iterative reconstruction acceleration method and system based on artificial intelligence | |
| CN105590332A (en) | Rapid algebraic reconstruction technique applied to computed tomography imaging | |
| CN117197349A (en) | CT image reconstruction method and device | |
| CN106707243A (en) | Generalized ROMP (Regularized Orthogonal Matching Pursuit) method for reconstructing radar signals | |
| KR20200058295A (en) | Method and Device of High Magnetic Field Magnetic Resonance Image Synthesis | |
| CN116664396A (en) | Quick and high-precision spine image stitching method | |
| Song et al. | Optimizing ADMM and Over-Relaxed ADMM Parameters for Linear Quadratic Problems | |
| CN110047048B (en) | Phase recovery improved algorithm based on MSE (mean square error) optimization |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200814 |