CN108765509B - Rapid image reconstruction method for linear imaging system - Google Patents
Rapid image reconstruction method for linear imaging system Download PDFInfo
- Publication number
- CN108765509B CN108765509B CN201810496663.5A CN201810496663A CN108765509B CN 108765509 B CN108765509 B CN 108765509B CN 201810496663 A CN201810496663 A CN 201810496663A CN 108765509 B CN108765509 B CN 108765509B
- Authority
- CN
- China
- Prior art keywords
- vector
- projection
- hyperplane
- prior information
- iteration
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 238000003384 imaging method Methods 0.000 title claims abstract description 36
- 239000013598 vector Substances 0.000 claims abstract description 107
- 238000005457 optimization Methods 0.000 claims abstract description 55
- 238000012804 iterative process Methods 0.000 claims abstract description 4
- 230000001133 acceleration Effects 0.000 claims description 39
- 230000008569 process Effects 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 claims description 8
- 238000013178 mathematical model Methods 0.000 claims description 6
- 238000005259 measurement Methods 0.000 claims description 5
- 125000004432 carbon atom Chemical group C* 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 abstract description 24
- 238000013170 computed tomography imaging Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 238000002591 computed tomography Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000009659 non-destructive testing Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
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 invention discloses a rapid image reconstruction method for a linear imaging system, wherein a plurality of times of accelerated adjustment and prior information optimization operations are added in each round of circular projection operation of the traditional algebraic reconstruction method, the distance from an iterative solution vector to an optimal solution vector is reduced in each time of accelerated adjustment operation, the convergence speed of the algorithm is obviously improved, the extra storage space and the calculated amount added by the algorithm are small and can be ignored, the iterative process does not depend on the iteration variable of the previous time, and the rapid image reconstruction of the linear imaging system under the condition of incomplete projection is easy to realize by directly combining the prior information.
Description
Technical Field
The invention belongs to the technical fields of biomedical imaging, nondestructive testing, image reconstruction and the like, and particularly relates to a rapid image reconstruction method for a linear imaging system.
Background
The image information is a branch with the largest information content and the most abundant content in the information engineering subject field, and the modern imaging system integrates the technologies of optics, electronics, computers, machinery and the like to convert the information in the objective world into the image information in various ways, thereby greatly enriching the human visual world and having wide application in numerous fields of medical diagnosis, nondestructive testing, detection, investigation and the like. At present, the new technology development of modern imaging systems is mainly characterized by the comprehensive application in the aspects of digitalization, multifunction, multidimensional, informatization and the like, and the imaging function and performance depend on the used imaging algorithm to a great extent. Therefore, imaging algorithms are key to imaging systems and are a hot spot in the research of imaging technology.
For many linear imaging systems, a linear model is often used to model the system mathematically, including solving a large-scale system of linear equations, and analytical and iterative methods are used as the imaging algorithms. Taking Computed Tomography (CT) imaging as an example, the current CT imaging system mostly uses a filtered back projection method (FBP), which has the advantages of small calculation amount, fast reconstruction speed, and good reconstructed image quality when the projection data is complete. However, in actual detection, it is often difficult to detect complete projection data for a number of objective reasons. When the projection data is insufficient, the quality of the image reconstructed by the FBP algorithm is seriously reduced, and the actual application requirements are difficult to meet.
Different from an analytic method represented by FBP, an Algebraic Reconstruction Technology (ART) converts an image reconstruction problem into a problem of solving a (large-scale) linear equation set at the initial stage of imaging, can select different objective functions to carry out iterative solution according to different requirements by combining prior information, and can reconstruct a high-quality image under the condition of sparse or incomplete projection data. However, the greatest disadvantage of the algebraic reconstruction method is that the calculated amount is large, the reconstruction speed is slow, and the time cost of the algorithm is an important factor which hinders the popularization of the algebraic reconstruction method.
In recent years, the progress of computer technology creates a new hardware environment for image reconstruction, and the technical advantages of the algebraic reconstruction method are more prominent. In order to enable an algebraic reconstruction algorithm to be widely applied to an actual imaging system, a simple and efficient fast reconstruction algorithm is urgently needed to be researched and developed, and fast and accurate imaging is realized under the condition of incomplete projection data by combining image prior information.
As the closest prior art, the patent "a fast algebraic reconstruction method applied to CT imaging" (application No. 201510988996.6) that the inventor has applied for discloses that on the basis of the traditional algebraic reconstruction method, partial hyperplanes are arbitrarily selected, and on these hyperplanes, the projection solution vectors obtained by the traditional algebraic reconstruction method are further subjected to accelerated adjustment operation, so that the solution vectors after adjustment and update are optimal on the connecting line of the original projection solution vectors and the corresponding vectors on the hyperplane in the previous iteration. The method obviously improves the convergence rate of the algorithm, obtains beneficial effects, but still has some defects when being directly applied to image reconstruction of incomplete projection data: the acceleration adjustment operation in the algorithm depends on the relationship between the current vector and the vector on the same hyperplane of the previous iteration, and if the prior information optimization operation is executed based on the current vector obtained by each iteration in combination with the prior information, the disturbance or change of the vector increased by the prior information optimization operation will destroy the dependency relationship, so that the subsequent iteration calculation of the algorithm is invalid. Therefore, the method disclosed in the applied patent (application No. 201510988996.6) is limited when combined with other a priori information optimization algorithms, and requires special processing, which is not favorable for the algorithm to flexibly combine a priori information to realize rapid image reconstruction in the case of incomplete projection data.
Disclosure of Invention
In order to solve the problems, the invention provides a rapid image reconstruction method for a linear imaging system, which directly performs multiple times of accelerated adjustment and prior information optimization operation in one iteration, only introduces one auxiliary vector and auxiliary scalar, hardly increases algorithm memory space and calculated amount additionally, does not depend on the last iteration variable in the iteration process, is easy to directly combine with prior information to realize rapid image reconstruction under the condition of incomplete data, introduces relaxation parameters and deduces corresponding formulas of auxiliary scalar updating and accelerated adjustment operation, and is suitable for noise environment. The technical scheme adopted by the invention is as follows:
a fast image reconstruction method for a linear imaging system includes the steps of firstly establishing a discretization mathematical model and a priori information mathematical model of the linear imaging system, converting an imaging problem of the linear imaging system into a large-scale linear equation set solving and priori information optimizing problem, then combining an acceleration adjusting operator and a priori information optimizing operator to iteratively solve the large-scale linear equation set, optimizing priori information, and reconstructing a high-quality image under the condition of sparse or incomplete projection data, wherein the solved large-scale linear equation set is represented by Ax b, x is vector representation of an image to be solved and is an N × 1 vector, two-dimensional or three-dimensional images are rearranged in one dimension, A is an M × N system matrix, and a vector a is used for rearranging a vector a to obtain the image to be reconstructed in one dimensioniRepresenting the transpose of the ith row vector of matrix A, b being the M × 1 measurement data vector, using biThe ith element representing vector b; it seeks toThe solution process comprises the following steps:
the method comprises the following steps: initialization: setting an initial solution vector x0Relaxation parameters rho and gamma, and the number T of operations for accelerating adjustment and prior information optimization in each loop iteration, and let the iteration number k be 0, and the current vector (i.e. the current estimation solution vector) x be x0Auxiliary scalar d is 0 and auxiliary vector z is x0；
Step two: and (3) an iterative process: repeating the following steps for the kth iteration until a convergence condition is satisfied:
(2.1) setting the order of the hyperplane projection by slA sequence number of the ith projection hyperplane is represented, wherein l is 1, 2.. M, M represents the number of hyperplanes, and t is 1, t represents the t-th acceleration adjustment and prior information optimization operation in the iteration;
(2.2) for each hyperplane i ═ slM, updating the current vector x by sequentially performing the following operations:
(2.2.1) performing a projection operation on the current vector x into the hyperplane i, updating the current vector x with the corresponding auxiliary scalar d:
[x,d]←proj(x,d,ρ,ai,bi)
wherein, [ x ', d']＝proj(x,d,ρ,ai,bi) Representing the operation of projection and updating an auxiliary scalar d, and using "←" as an assignment symbol;
(2.2.2) if
Then, further performing acceleration adjustment and prior information optimization operations on the current vector x, updating the current vector x, and the corresponding auxiliary vector z and auxiliary scalar d:
[x,z,d]←lineAcc(x,z,d,γ)
t←t+1；
in the above formula, the first and second carbon atoms are,
represents rounding up; [ x ', z ', d ']lineAcc (x, z, d, γ) represents an acceleration adjustment and a priori information optimization operation;
(2.3)k←k+1。
furthermore, in the step one, the value ranges of the relaxation parameters ρ and γ are (0, 1), and the value range of the accelerated adjustment and prior information optimization operation number T is more than or equal to 1 and less than or equal to T and less than or equal to M.
Further, in the step (2.1), a hyperplane projection sequence is set by adopting a fixed projection sequence or a random projection sequence.
Further, in the step (2.1), if a hyperplane projection order is set by using a fixed projection order, the number i ═ s of the l (i ═ 1, 2.. M) th projection hyperplane is iterated each timelIf the value of (2) is fixed, the acceleration adjustment and the prior information optimization operation in each iteration step (2.2.2) aim at the same hyperplane serial number; on the contrary, if the projection sequence of the hyperplane is set by adopting the random projection sequence, the serial number i of the first projection hyperplane in each iteration is equal to slIs random, the acceleration adjustment and the prior information optimization operation in each iteration step (2.2.2) are for different hyperplane numbers.
Further, the projection operation [ x ', d ' in the step (2.2.1) ']＝proj(x,d,ρ,ai,bi) The definition is as follows:
[x',d']＝proj(x,d,ρ,ai,bi)
λ＝bi-<ai,x>(1)
when ρ ═ 1, projection operation [ x ', d']＝proj(x,d,λ,ai,bi) Is the projected point of the current vector x to the hyperplane i. After the projection operation is executed, the following requirements are met: | x*-x||2-||x*-x'||2D-d', where x*Is the true solution of the equation set, i.e., the intersection of all hyperplanes.
Further, the acceleration adjustment and prior information optimization operation [ x ', z ', d ' ] — lineAcc (x, z, d, γ) in step (2.2.2) is defined as follows:
[x',z',d']＝lineAcc(x,z,d,γ)
u＝x-z (4)
x″＝x+γβu (6)
x'＝Sup(x″) (7)
z'＝x' (8)
d'＝0 (9)
acceleration adjustment and prior information optimization operation [ x ', z ', d ']The equations (4), (5), (6) of lineAcc (x, z, d, γ) perform the acceleration adjustment operation when | | x*-z||2-||x*-x||2When the formula is 0-d (which is satisfied until the algorithm iteration process performs the acceleration adjustment and prior information optimization operation) and γ is 1, x "obtained by the formula (6) is x" which is a real solution x to the equation system on a connecting line between the vector x and the corresponding auxiliary vector z*The closest point; when gamma is more than 0 and less than or equal to 1, x' is more than the real solution x of the vector x-ray equation set before updating*Closer. Wherein x is*Is the intersection of all hyperplanes, i.e. for i ═ 1,2, …, M, there are:<ai,x*>＝bi。[x',z',d']equation (7) of lineAcc (x, z, d, γ) performs a priori information optimization operation, and x' returned by x ═ Sup (x ″) is closer to the priori information of the image than x ″. There are many representations of the prior information of an image, many times described in some mathematical transformation. For example, an effective and widely used mathematical description is the Total Variation (TV) norm of an image, which is relatively small. When the TV norm is used to describe the prior information, the prior information optimization operation x '═ Sup (x ″) in the method of the present invention is to reduce the TV norm of the image, i.e. x' is returned with a smaller TV norm of the image than x ″. [ x ', z ', d ']Equations (8) and (9) of lineAcc (x, z, d, γ) update the auxiliary vector z and the auxiliary scalar d, respectively. Performing accelerated adjustment and prior information optimization operation [ x ', z ', d ']After lineAcc (x, z, d, γ), in general, the updated vector x' is larger than the vector before updatex has less measurement error and is closer to the prior information of the image (e.g., the TV norm is smaller).
The theoretical analysis of the invention is as follows:
the basic idea of algebraic reconstruction techniques is: considering each equation as a hyperplane of an N-dimensional space, from an initial solution vector x0Starting from the above, in each iteration, according to a certain projection sequence, sequentially performing one round of circular projection operation on each hyperplane to gradually approach to the real solution of the equation set. Taking sequential projection (sequential projection hyperplane numbers are 1,2, …, M), starting from a current estimation solution vector (current vector for short) x, at each iteration, the following operations are sequentially performed according to the hyperplane projection sequence i being 1,2, …, M, and the estimation solution vector x is continuously updated:
λ←bi-<ai,x>(10)
in the above formula, < > is the inner product symbol, and "←" is the assignment symbol. Equation (11) describes the projection operation, where the vector x on the right side of equation is the vector before the projection operation is performed, and the vector x on the left side of equation is the vector updated after the projection operation is performed. ρ is a relaxation parameter of the ART algorithm, and is applied to a noise environment, and a value range is generally between (0, 1), when ρ is equal to 1, equations (10) and (11) project a current vector x onto a hyperplane (abbreviated as hyperplane i) with a sequence number i.
In the method, if the acceleration adjustment and the prior information optimization operation in the step (2.2.2) are removed, or the number T of the acceleration adjustment and the prior information optimization operation is set to be 0, the method is degenerated to the traditional algebraic reconstruction method. After the acceleration adjustment and the prior information optimization operation in the step (2.2.2) are added, the algorithm speed is improved, and meanwhile, a high-quality image is reconstructed under the condition of sparse or incomplete projection data.
First, it is explained that in step (2.2.1), the projection operation [ x ', d ' to the hyperplane i is performed on the current vector x ']＝proj(x,d,ρ,ai,bi) After that, there are:
||x*-x||2-||x*-x'||2＝d-d' (12)
the following was demonstrated:
when the relaxation parameter ρ in the formula (2) is 1, a vector x' (expressed as x) obtained by the formula (2) is obtainedp) For the projection point of the current vector x on the hyperplane i:
i.e. xpSatisfies the following conditions:
<x-xp,x*-xp>＝0 (14)
the method comprises the following steps:
||x*-x||2＝||x*-xp||2+||x-xp||2(15)
||x*-x'||2＝||x*-xp||2+||x'-xp||2(16)
the following equations (2) and (13) can be obtained:
the bond formulas (15), (16), (17), (18) give:
the following equations (19) and (3) can be obtained:
||x*-x||2-||x*-x'||2＝d-d'
and (5) finishing the certification.
During the iteration of the algorithm, vector x corresponds to an auxiliary scalar d and auxiliary vector z corresponds to an auxiliary scalar 0. In combination with the above properties, it can be seen that the following condition is satisfied before the acceleration adjustment and prior information optimization operations are performed in step (2.2.2):
||x*-z||2-||x*-x||2＝0-d (20)
it is explained that if γ ═ 1, the acceleration adjustment and prior information optimization operation [ x ', z', d 'is performed in step (2.2.2)']After equation (6) of lineAcc (x, z, d, γ), the returned x "is the true solution x to the system of equations on the connecting line of the vector x and the corresponding auxiliary vector z*The closest point, i.e. satisfy
<u，x*-x″>＝0 (21)
In the above equation, u and x "are defined by equations (4) and (6), respectively, (where γ is equal to 1), that is, if x" returned by equation (6) for updating x is x "which is a real solution x from the equation set on a connecting line between the vector x and the corresponding auxiliary vector z when γ is equal to 1*In the most recent case, β in equation (6) is defined by equation (5).
The following was demonstrated:
the following are obtained by combining the formulas (21), (4) and (6):
||x*-x||2＝||x*-x″||2+β2||u||2(22)
||x*-z||2＝||x*-x″||2+(1+β)2||u||2(23)
from the formulae (22) and (23), it is understood that:
||x*-z||2-||x*-x||2＝(2β+1)||u||2(24)
the following are shown by the combined formulae (20) and (24):
namely, the calculation expression of the formula (5).
And (5) finishing the certification.
In consideration of noise, relaxation parameters ρ, γ are introduced into equations (2) and (6). When gamma is more than 0 and less than or equal to 1, optimizing operation [ x ', z', d 'by acceleration adjustment and prior information']The returned vector x "of equation (6) is greater than the vector x before updating from the true solution x of the system of equations*Closer.
The beneficial effects of the technical scheme are as follows:
1. compared with the prior art (namely the method disclosed by the inventor's previous patent application (application number 201510988996.6)), the method disclosed by the invention directly performs multiple times of acceleration adjustment and prior information optimization operation within one loop iteration, and because the iteration process does not depend on the previous iteration variable, the rapid image reconstruction under the incomplete projection condition can be realized more easily by directly combining the prior information.
2. Compared with the method disclosed by the prior art (namely the method disclosed by the previous patent application (application number 201510988996.6) of the inventor), the method appropriately introduces relaxation parameters to the noise environment and derives corresponding formulas of auxiliary scalar updating and accelerating adjustment operation, so that the application range of the method is expanded.
3. Improving the convergence speed of an image reconstruction algorithm: compared with the traditional algebraic reconstruction method, the method only adds one auxiliary vector z, and the added extra storage space can be completely ignored. In practical application, the number T of the acceleration adjustment and prior information optimization operations in each iteration can be set to be small (for example, less than 10), so that the extra calculation amount can be ignored, and the simulation result shows that the convergence rate of the algorithm can be significantly improved.
4. The imaging system application can be effectively expanded: the method can be applied to different linear imaging systems, flexibly utilizes corresponding prior information, and realizes quick image reconstruction under the condition of incomplete data.
Drawings
The present invention will be described in further detail with reference to the accompanying drawings.
FIG. 1 is a flow chart of a fast image reconstruction method for a linear imaging system according to the present invention;
FIG. 2 is a diagram of the method of the present invention for a hyperplane i ═ slA flow diagram illustrating the update operation performed on the current vector x;
FIG. 3 is a diagram of a comparison of a reconstructed image and an original image according to an embodiment of the method of the present invention;
fig. 4 is a graph of the convergence of square vector error for an embodiment of the method of the present invention and an algorithm incorporating ART and TV norm optimization.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be further described with reference to the accompanying drawings and specific embodiments.
In this embodiment, referring to fig. 1 and 2, a fast image reconstruction method for a CT imaging system includes first establishing a discretization mathematical model and a priori information mathematical model of the CT imaging system, converting the CT imaging problem into a large-scale linear equation set solving and a priori information optimizing problem, then iteratively solving the large-scale linear equation set by combining an accelerated modification operator and a priori information optimizing operator, optimizing priori information, reconstructing a high-quality image under a sparse or incomplete projection data condition, where the solved large-scale linear equation set is represented by Ax ═ b, where x is a vector representation of an image to be solved and is an N × 1 vector, and rearranging a two-dimensional images in one dimension, where a is an M × N system matrix and is a vector aiRepresenting the transpose of the ith row vector of matrix A, b being the M × 1 measurement data vector, using biThe ith element representing vector b; the solving process comprises the following steps:
the method comprises the following steps: initialization: setting an initial solution vector x0Relaxation parameters rho and gamma, and the number T of operations for accelerating adjustment and prior information optimization inside each loop iteration, and the number k of iterations is 0, and the current vector x is x0Auxiliary scalar d is 0 and auxiliary vector z is x0；
In the first step, the value ranges of the relaxation parameters rho and gamma are (0, 1), and the value range of the accelerated adjustment and prior information optimization operation number T is more than or equal to 1 and less than or equal to M.
Step two: and (3) an iterative process: repeating the following steps for the kth iteration until a convergence condition is satisfied:
(2.1) setting the order of the hyperplane projection by slA sequence number of the ith projection hyperplane is represented, wherein l is 1, 2.. M, M represents the number of hyperplanes, and t is 1, t is used for representing the t-th acceleration adjustment and prior information optimization operation in the iteration;
and (3) setting a hyperplane projection sequence by adopting a fixed projection sequence or a random projection sequence in the step (2.1).
(2.2) for each hyperplane i ═ slM, updating the current vector x by sequentially performing the following operations:
(2.2.1) performing a projection operation on the current vector x into the hyperplane i, updating the current vector x with the corresponding auxiliary scalar d:
[x,d]←proj(x,d,ρ,ai,bi)
wherein, [ x ', d']＝proj(x,d,ρ,ai,bi) Representing the operation of projection and updating an auxiliary scalar d, and using "←" as an assignment symbol;
the projection operation [ x ', d ' in the step (2.2.1) ']＝proj(x,d,ρ,ai,bi) The definition is as follows:
[x',d']＝proj(x,d,ρ,ai,bi)
λ＝bi-<ai,x>(1)
when ρ ═ 1, projection operation [ x ', d']＝proj(x,d,λ,ai,bi) Is the projected point of the current vector x to the hyperplane i. After the projection operation is executed, the following requirements are met: | x*-x||2-||x*-x'||2D-d', where x*Is a true solution of the system of equations.
(2.2.2) if
Then, further performing acceleration adjustment and prior information optimization operations on the current vector x, updating the current vector x, and the corresponding auxiliary vector z and auxiliary scalar d:
[x,z,d]←lineAcc(x,z,d,γ)
t←t+1；
in the above formula, the first and second carbon atoms are,
represents rounding up; [ x ', z ', d ']lineAcc (x, z, d, γ) represents an acceleration adjustment and a priori information optimization operation;
the acceleration adjustment and prior information optimization operation [ x ', z ', d ' ] — lineAcc (x, z, d, γ) in step (2.2.2) is defined as follows:
[x',z',d']＝lineAcc(x,z,d,γ)
u＝x-z (4)
x″＝x+γβu (6)
x'＝Sup(x″) (7)
z'＝x' (8)
d'＝0 (9)
acceleration adjustment and prior information optimization operation [ x ', z ', d ']The equations (4), (5), (6) of lineAcc (x, z, d, γ) perform the acceleration adjustment operation when | | x*-z||2-||x*-x||2When the formula is 0-d (which is satisfied until the algorithm iteration process performs the acceleration adjustment and prior information optimization operation) and γ is 1, x "obtained by the formula (6) is x" which is a real solution x to the equation system on a connecting line between the vector x and the corresponding auxiliary vector z*The closest point; when gamma is more than 0 and less than or equal to 1, x' is more than the real solution x of the vector x-ray equation set before updating*Closer. [ x ', z ', d ']Equation (7) of lineAcc (x, z, d, γ) performs a priori information optimization operation, and x' returned by x ═ Sup (x ") is generally closer to the priori information of the image than x ″. There are many representations of the prior information of an image, and many times it is described by some mathematical transformation, such as a Total Variation (TV) norm. When the TV norm is used to describe the prior information, the prior information optimization operation x '═ Sup (x ″) in the method of the present invention is to reduce the TV norm of the image, i.e. x' is returned with a smaller TV norm of the image than x ″. [ x ', z ', d ']Equations (8) and (9) of lineAcc (x, z, d, γ) update the auxiliary vector z and the auxiliary scalar d, respectively. Performing accelerated adjustment and prior information optimization operation [ x ', z ', d ']After lineAcc (x, z, d, γ), in general termsThe updated vector x' has a smaller measurement error than the vector x before updating, and is closer to the prior information of the image (e.g., the TV norm is smaller).
In the step (2.1), a fixed projection sequence is adopted to set a projection sequence of the hyperplane, and the sequence number i of the ith (l is 1, 2.. M) projection hyperplane is s in each iterationlIf the value of (2) is fixed, the acceleration adjustment and the prior information optimization operation in the next step (2.2.2) in each iteration are directed to the same hyperplane serial number; on the contrary, if the projection sequence of the hyperplane is set by adopting the random projection sequence, the serial number i of the first projection hyperplane in each iteration is equal to slIs random, each iteration the acceleration adjustment and prior information optimization operations in the following step (2.2.2) are for different hyperplane numbers.
(2.3) k ← k +1, i.e. the number of iterations plus 1, ready for the next iteration.
The invention is further illustrated below with reference to 1 specific example.
In this embodiment, corresponding to the situation that the projection data is not complete, the CT image to be reconstructed is a 128 × 128 two-dimensional image, i.e. the original image in fig. 3, corresponding to the vector x represented as 16384 × 1*A is the system matrix with the dimension 7200 × 16384, i.e., M7200, N16384, b Ax*+ e is the measured data vector of 16384 × 1 (matrix a contains a small number of all-zero rows that are rejected during processing). e is a zero-mean gaussian white noise vector with an amplitude intensity of 1% of measured data b.
Under the incomplete projection condition, in order to reconstruct a high-quality image, some a priori knowledge needs to be introduced for regularization. In this example, a Total Variation (TV) regularization technique is employed, and the problem of CT image reconstruction after regularization is converted into the minimization problem of the following objective function:
in this embodiment, η is 0.004, the parameters of the method of the present invention are set as ρ 1, γ 1, and the initial vector x in the initialization phase0Set to all zero vectors, each iteration of the loopAnd (4) carrying out the operation number T of acceleration adjustment and prior information optimization in the generation, and adopting a random projection sequence.
In each iteration process, firstly, a group of random permutation and combination of 1, 2.. M is generated as a projection sequence, and then projection and acceleration adjustment and prior information optimization operation are carried out on each hyperplane according to the randomly generated projection sequence.
For items 1800, 3600, 5400, 7200 (i.e., items)
) The other projection hyperplane only executes projection operation (2.2.1) and does not execute acceleration adjustment and prior information optimization operation (2.2.2); for the 1800 th, 3600 th, 5400 th and 7200 th projection hyperplanes, the acceleration adjustment and prior information optimization operation (2.2.2) is also performed after the projection operation (2.2.1) is performed. Because a random projection sequence is adopted, the serial numbers of the 1800 th, 3600 th, 5400 th and 7200 th projection hyperplanes in each iteration are generally different. For example, the number of the 1800 th projection hyperplane in the first iteration is 353, and the number of the 1800 th projection hyperplane in the second iteration is 3478, which means that the accelerated adjustment and prior information optimization operation is performed on the hyperplane with the number of 353 in the first iteration, and the accelerated adjustment and prior information optimization operation is performed on the hyperplane with the number of 3478 in the second iteration.
Because the method does not depend on the previous iteration variable in the iteration process, the rapid image reconstruction under the incomplete projection condition is easy to realize by directly combining the prior information. The comparison graph of the image reconstructed by the method of the invention and the original image is shown in fig. 3, the square vector error convergence curve graph of the algorithm combining ART and TV norm optimization by the method of the invention is shown in fig. 4, and the graph shows that the convergence speed is obviously improved by the method of the invention.
The foregoing shows and describes the general principles, essential features, and advantages of the invention. It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.
Claims (6)
1. A fast image reconstruction method for a linear imaging system is characterized by firstly establishing a discretization mathematical model and a priori information mathematical model of the linear imaging system, converting an imaging problem of the linear imaging system into a large-scale linear equation set solving and a priori information optimizing problem, then combining an accelerated adjustment operator and a priori information optimizing operator to iteratively solve the large-scale linear equation set, optimizing priori information, and reconstructing a high-quality image under the condition of sparse or incomplete projection data, wherein the solved large-scale linear equation set is represented by Ax b, x is vector representation of an image to be solved and is N × 1 vector, the image to be solved is obtained by rearranging two-dimensional or three-dimensional images according to one dimension, A is an M × N system matrix, and vector a is used for reconstructing the image to be solved according to one dimensioniRepresenting the transpose of the ith row vector of matrix A, b being the M × 1 measurement data vector, using biThe ith element representing vector b; the solving process comprises the following steps:
the method comprises the following steps: initialization: setting an initial solution vector x0Relaxation parameters rho and gamma, and the number T of operations for accelerating adjustment and prior information optimization inside each iteration, and the number k of iterations is 0, and the current vector x is x0Auxiliary scalar d is 0 and auxiliary vector z is x0；
Step two: and (3) an iterative process: repeating the following steps for the kth iteration until a convergence condition is satisfied:
(2.1) setting the order of the hyperplane projection by slA sequence number of the ith projection hyperplane is represented, wherein l is 1, 2.. M, M represents the number of hyperplanes, and t is 1, t represents the t-th acceleration adjustment and prior information optimization operation in the iteration;
(2.2) for each hyperplane i ═ slM, updating the current vector x by sequentially performing the following operations:
(2.2.1) performing a projection operation on the current vector x into the hyperplane i, updating the current vector x with the corresponding auxiliary scalar d:
[x,d]←proj(x,d,ρ,ai,bi)
wherein, [ x ', d']＝proj(x,d,ρ,ai,bi) Representing the updating operation of projection and an auxiliary scalar d, wherein x 'and d' in the formula are respectively the updating result of the current vector x and the auxiliary scalar d after the projection operation is executed, and "←" is an assignment symbol;
(2.2.2) if
Then, further performing acceleration adjustment and prior information optimization operations on the current vector x, updating the current vector x, and the corresponding auxiliary vector z and auxiliary scalar d:
[x,z,d]←lineAcc(x,z,d,γ)
t←t+1；
in the above formula, the first and second carbon atoms are,
represents rounding up; [ x ', z ', d ']Representing acceleration adjustment and prior information optimization operations, wherein x ', z ' and d ' in the formula are update results of a current vector x, an auxiliary vector z and an auxiliary scalar d after the acceleration adjustment and prior information optimization operations are performed respectively;
(2.3)k←k+1。
2. the fast image reconstruction method for a linear imaging system according to claim 1, wherein in the first step, the range of values of the relaxation parameters p, γ is (0, 1), and the range of values of the number T of the accelerated adjustment and prior information optimization operations is 1 ≤ T ≤ M.
3. The fast image reconstruction method for linear imaging system according to claim 1, characterized in that the step (2.1) sets the order of the hyperplane projections using a fixed projection order or a random projection order.
4. A method for linear imaging according to claim 3The fast image reconstruction method of the system is characterized in that, if the projection order of the hyperplane is set by adopting a fixed projection order in the step (2.1), the serial number i of the ith projection hyperplane is s in each iterationlIf the value of (2) is fixed, the acceleration adjustment and the prior information optimization operation in each iteration step (2.2.2) aim at the same hyperplane serial number; on the contrary, if the projection sequence of the hyperplane is set by adopting the random projection sequence, the serial number i of the first projection hyperplane in each iteration is equal to slIs random, the acceleration adjustment and the prior information optimization operation in each iteration step (2.2.2) are for different hyperplane numbers.
5. A fast image reconstruction method for linear imaging system according to claim 1, characterized in that the projection and auxiliary scalar d update operation [ x ', d ' in said step (2.2.1) ']＝proj(x,d,ρ,ai,bi) The definition is as follows:
[x',d']＝proj(x,d,ρ,ai,bi)
in the above formula, < - > represents the inner product operation of the vector.
6. A fast image reconstruction method for a linear imaging system according to claim 1, characterized in that the acceleration adjustment and prior information optimization operation [ x ', z ', d ' ] -lineAcc (x, z, d, γ) in step (2.2.2) is defined as follows:
[x',z',d']＝lineAcc(x,z,d,γ)
x′＝Sup(x″)
z'＝x'
d'＝0；
wherein x' is the real solution x of the equation set on the connecting line of the vector x and the corresponding auxiliary vector z*The closest point.
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 CN108765509A (en) | 2018-11-06 |
| CN108765509B true 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) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115619890B (en) * | 2022-12-05 | 2023-04-07 | 山东省计算中心（国家超级计算济南中心） | Tomography method and system for solving linear equation set based on parallel random iteration |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8055095B2 (en) * | 2008-01-23 | 2011-11-08 | Sparsense, Inc. | Parallel and adaptive signal processing |
| US9031297B2 (en) * | 2013-02-01 | 2015-05-12 | Kabushiki Kaisha Toshiba | Alternative noise map estimation methods for CT images |
-
2018
- 2018-05-22 CN CN201810496663.5A patent/CN108765509B/en not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
|---|
| Algebraic Reconstruction Technique with adaptive relaxation parameter based on hyperplane distance and data noise level;Chuan Lin 等;《2016 IEEE MTT-S International Conference on Numerical Electromagnetic and Multiphysics Modeling and Optimization (NEMO)》;20160729;全文 * |
| 基于不完备投影数据重建的四种迭代算法比较研究;阙介民 等;《CT理论与应用研究》;20120615;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108765509A (en) | 2018-11-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sun et al. | Block coordinate regularization by denoising | |
| CN110047113B (en) | Neural network training method and apparatus, image processing method and apparatus, and storage medium | |
| Sreehari et al. | Plug-and-play priors for bright field electron tomography and sparse interpolation | |
| Zhuge et al. | TVR-DART: A more robust algorithm for discrete tomography from limited projection data with automated gray value estimation | |
| Chen et al. | Adaptively regularized constrained total least-squares image restoration | |
| CN108765514B (en) | Acceleration method and device for CT image reconstruction | |
| Zhang et al. | JSR-Net: a deep network for joint spatial-radon domain CT reconstruction from incomplete data | |
| Li et al. | Reconstruction of heterogeneous materials via stochastic optimization of limited-angle X-ray tomographic projections | |
| CN111899314A (en) | Robust CBCT reconstruction method based on low-rank tensor decomposition and total variation regularization | |
| Li et al. | Sparse CT reconstruction based on multi-direction anisotropic total variation (MDATV) | |
| Kadu et al. | A convex formulation for binary tomography | |
| CN108765509B (en) | Rapid image reconstruction method for linear imaging system | |
| Wang et al. | Mars image super-resolution based on generative adversarial network | |
| Zhang et al. | Nonconvex log-sum function-based majorization–minimization framework for seismic data reconstruction | |
| CN105608719B (en) | A kind of rapid CT image rebuilding method based on two benches projection adjustment | |
| Hashemi et al. | Fast fan/parallel beam CS-based low-dose CT reconstruction | |
| CN117408991A (en) | Image anomaly detection method and device based on reconstruction resistance and storage medium | |
| Han et al. | Weighting algorithm and relaxation strategies of the landweber method for image reconstruction | |
| US8379948B2 (en) | Methods and systems for fast iterative reconstruction using separable system models | |
| CN114913262B (en) | Nuclear magnetic resonance imaging method and system with combined optimization of sampling mode and reconstruction algorithm | |
| Kudo et al. | Metal artifact reduction in CT using fault-tolerant image reconstruction | |
| CN110060314B (en) | CT iterative reconstruction acceleration method and system based on artificial intelligence | |
| Liu et al. | Accelerated augmented Lagrangian method for total variation minimization | |
| Barrett | Perspectives on SPECT | |
| Di Serafino et al. | TGV-based restoration of Poissonian images with automatic estimation of the regularization parameter |
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 |