CN108765509B - Rapid image reconstruction method for linear imaging system - Google Patents

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

A Fast Image Reconstruction Method for Linear Imaging Systems

technical field

The invention belongs to the technical fields of biomedical imaging, non-destructive testing, image reconstruction and the like, in particular to a fast image reconstruction method for a linear imaging system.

Background technique

Image information is a branch with the largest amount of information and the most abundant content in the field of information engineering. The modern imaging system integrates optics, electronics, computers, machinery and other technologies to transform the information in the objective world into image information in various ways. It has enriched the human visual world and has been widely used in many fields such as medical diagnosis, non-destructive testing, detection, and reconnaissance. At present, the characteristics of the development of new technologies in modern imaging systems are mainly reflected in the comprehensive application of digitization, multi-functionalization, multi-dimensionalization and informatization. The functions and performance of imaging largely depend on the imaging algorithms used. Therefore, the imaging algorithm is the key to the imaging system, and it is also the hotspot of imaging technology research.

For many linear imaging systems, a linear model can often be used to model mathematically, which includes the solution of large-scale linear equations. The imaging algorithms that can be used include analytical and iterative solutions. Taking computed tomography (CT) imaging as an example, the current CT imaging systems mostly use the filtered back projection (FBP) method, 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, due to many objective reasons, it is often difficult to detect complete projection data. When the projection data is insufficient, the quality of the image reconstructed by the FBP algorithm is seriously degraded, and it is difficult to meet the needs of practical applications.

Different from the analytical method represented by FBP, Algebraic Reconstruction Technology (ART) transforms the image reconstruction problem into the problem of solving a (large-scale) linear equation system at the beginning of imaging, which can be combined with prior information to select different targets according to different requirements. The function is iteratively solved, and high-quality images can still be reconstructed under the condition of sparse or incomplete projection data. But the biggest disadvantage of the algebraic reconstruction method is the large amount of calculation and the slow reconstruction speed. The time cost of the algorithm is an important factor hindering its popularization.

In recent years, the advancement of computer technology has created a new hardware environment for image reconstruction, and the technical advantages of algebraic reconstruction methods have become more prominent. In order to make the algebraic reconstruction algorithm widely used in the actual imaging system, it is urgent to research and develop a simple and efficient fast reconstruction algorithm, combining with the prior information of the image, to achieve fast and accurate imaging in the case of incomplete projection data.

As the closest prior art, the patent "A Fast Algebraic Reconstruction Method Applied to CT Imaging" (application number 201510988996.6) that the inventor has applied for discloses that on the basis of the traditional algebraic reconstruction method, some hyperplanes are arbitrarily selected, On these hyperplanes, the projection solution vector obtained by the traditional algebraic reconstruction method is further accelerated and adjusted, so that the adjusted and updated solution vector is optimal on the connecting line between the original projection solution vector and the corresponding vector on the hyperplane in the previous iteration. This method significantly improves the convergence speed of the algorithm and obtains beneficial effects, but there are still some shortcomings when it is directly applied to image reconstruction of incomplete projection data: the accelerated adjustment operation in the algorithm depends on the current vector on the same hyperplane as the previous iteration. If the prior information optimization operation is performed on the basis of the current vector obtained by each iteration and the prior information is combined, the disturbance or change to the vector due to the prior information optimization operation will destroy this dependency. The subsequent iterative calculation of the algorithm fails. Therefore, the method disclosed in the applied patent (application number 201510988996.6) is subject to certain limitations when combined with other prior information optimization operators, and requires special processing, which is not conducive to the algorithm flexibly combining prior information Enables fast image reconstruction with incomplete projection data.

SUMMARY OF THE INVENTION

In order to solve the above problems, the present invention provides a fast image reconstruction method for a linear imaging system. The method directly performs multiple acceleration adjustments and prior information optimization operations within one iteration, and only introduces an auxiliary vector and auxiliary scalar, which is almost It does not increase the storage and calculation amount of the algorithm, the iteration process does not depend on the previous iteration variables, and it is easy to directly combine the prior information to achieve fast image reconstruction in the case of incomplete data. The method introduces relaxation parameters and derives the auxiliary scalar update and Corresponding formula for speeding up adjustment operations, suitable for noisy environments. The technical scheme adopted in the present invention is:

A fast image reconstruction method for linear imaging systems. First, a discrete mathematical model and a priori information mathematical model of the linear imaging system are established, and the imaging problem of the linear imaging system is transformed into a large-scale linear equation system solution and prior information optimization. Then combine the accelerated adjustment operator and the prior information optimization operator to iteratively solve the large-scale linear equation system, optimize the prior information, and reconstruct high-quality images under the condition of sparse or incomplete projection data; the solved large-scale linear equation system It is expressed as Ax=b, where x is the vector representation of the image to be obtained, which is an N×1 vector, obtained by rearranging a two-dimensional or three-dimensional image in one dimension, A is an M×N system matrix, represented by a vector a _i The transposition of the vector in the ith row of matrix A, b is the M×1 measurement data vector, and the ith element of vector b is represented by b _i ; the solution process includes the following steps:

Step 1: Initialization: Set the initial solution vector x ⁰ , relaxation parameters ρ, γ, and the number T of acceleration adjustment and prior information optimization operations performed inside each loop iteration, and set the number of iterations k=0, the current vector ( That is, the current estimated solution vector) x=x ⁰ , the auxiliary scalar d=0, and the auxiliary vector z=x ⁰ ;

Step 2: Iterative process: Repeat the following steps for the k-th iteration until the convergence conditions are met:

(2.1) Set the hyperplane projection order, and use s _l to denote the serial number of the l-th projection hyperplane, where l=1, 2,...M, M denotes the number of hyperplanes, and let t=1, t represents the t-th time to perform acceleration adjustment and prior information optimization operations within the iteration;

(2.2) For each hyperplane i=s _l , l=1, 2,...M, perform the following operations to update the current vector x in sequence:

(2.2.1) Perform the projection operation on the current vector x to the hyperplane i, and update the current vector x and the corresponding auxiliary scalar d:

[x,d]←proj(x,d,ρ,a _i ,b _i )

Among them, [x', d']=proj(x, d, ρ, a _i , b _i ), represents the update operation of projection and auxiliary scalar d, and "←" is the assignment symbol;

则对当前向量x进一步执行加速调整与先验信息优化操作，更新当前向量x，及对应的辅助向量z与辅助标量d：(2.2.2) If

Then further perform acceleration adjustment and prior information optimization operations on the current vector x, update the current vector x, and the corresponding auxiliary vector z and auxiliary scalar d:

[x,z,d]←lineAcc(x,z,d,γ)

t←t+1;

表示向上取整；[x',z',d']＝lineAcc(x,z,d,γ)表示加速调整与先验信息优化操作；In the above formula,

Represents rounding up; [x', z', d']=lineAcc(x, z, d, γ) represents accelerated adjustment and prior information optimization operations;

(2.3) k←k+1.

Further, in the first step, the value range of the relaxation parameters ρ and γ is (0, 1], and the value range of the number T of acceleration adjustment and prior information optimization operations is 1≤T≤M.

Further, in the step (2.1), a fixed projection order or a random projection order is used to set the hyperplane projection order.

Further, in the step (2.1), if a fixed projection sequence is used to set the hyperplane projection sequence, the sequence number i=s _l of the lth (l=1, 2,...M) projection hyperplane in each iteration is fixed, the acceleration adjustment and prior information optimization operation in each iteration step (2.2.2) are aimed at the same hyperplane number; otherwise, if the random projection order is used to set the hyperplane projection order, each The sequence number i=s _l of the l-th projection hyperplane in the second iteration is random, and the acceleration adjustment and prior information optimization operations in each iteration step (2.2.2) are aimed at different hyperplane numbers.

Further, the projection operation [x', d']=proj(x, d, ρ, a _i , b _i ) in the step (2.2.1) is defined as follows:

[x',d']=proj(x,d,ρ,a _i ,b _i )

λ=b _i -<a _i ,x> (1)

When ρ=1, the output x' of the projection operation [x', d']=proj(x, d, λ, a _i , b _i ) is the projection point of the current vector x to the hyperplane i. After performing the projection operation, satisfy: ||x ^* -x|| ² -||x ^* -x'|| ² =d-d', where x ^* is the true solution of the system of equations, that is, the intersection of all hyperplanes .

Further, the acceleration adjustment and prior information optimization operation [x', z', d']=lineAcc(x, z, d, γ) in the 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)

Accelerated adjustment and prior information optimization operation [x', z', d']=lineAcc(x, z, d, γ) Formulas (4), (5), (6) perform the accelerated adjustment operation, when || x ^* -z|| ² -||x ^* -x|| ² =0-d (this formula is always satisfied before the acceleration adjustment and prior information optimization operation is performed in the iterative process of the algorithm) and when γ=1, formula (6 ) obtained x″ is the closest point to the real solution x ^* of the system of equations on the connection line between the vector x and the corresponding auxiliary vector z; when 0<γ≤1, x″ is closer than the vector x before the update to the real solution x of the system of equations ^* is closer. where x ^* is the intersection of all hyperplanes, ie for i=1,2,...,M, there are: <a _i ^,x*> =b _i . [x', z', d']=lineAcc(x, z, d, γ) formula (7) performs prior information optimization operation, x'=Sup(x") returns x' ratio x" and the image The prior information is closer. The prior information of the image has many representations, and is often described by some mathematical transformation. For example, an effective and widely used mathematical description method is the Total Variation (TV) norm of the image, and the TV norm of the general image 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, that is, the returned x' has a smaller value than x" Image TV norm. Formulas (8) and (9) of [x', z', d']=lineAcc(x, z, d, γ) update the auxiliary vector z and the auxiliary scalar d, respectively. After performing the accelerated adjustment and prior information optimization operation [x', z', d'] = lineAcc(x, z, d, γ), in general, the updated vector x' has more Small measurement error, which is closer to the prior information of the image (eg, the TV norm is smaller).

The theoretical analysis of the present invention is:

The basic idea of the algebraic reconstruction technique is to regard each equation as a hyperplane in an N-dimensional space, starting from the initial solution vector x ⁰ , and in each iteration, according to a certain projection sequence, perform a round of cyclic projection to each hyperplane in turn. operation, step-by-step approximation to the true solution of the system of equations. Taking sequential projection (the hyperplane numbers of the sequential projection are 1, 2, ..., M) as an example, in each iteration, starting from the current estimated solution vector (referred to as the current vector) x, according to the hyperplane projection order i=1, 2 ,...,M, perform the following operations in sequence to continuously update the estimated solution vector x:

λ←b _i -<a _i ,x> (10)

<·> in the above formula is the inner product symbol, and "←" is the assignment symbol. Equation (11) describes the projection operation, the vector x on the right side of the expression is the vector before the projection operation is performed, and the vector x on the left side of the expression is the vector updated after the projection operation is performed. ρ is the relaxation parameter of the ART algorithm, which is applied to the noise environment. The value range is generally between (0, 1]. When ρ = 1, formulas (10) and (11) project the current vector x to the sequence number i. on the hyperplane (referred to as hyperplane i).

In the method of the present invention, if the accelerated adjustment and prior information optimization operations in step (2.2.2) are removed, or the number T of the accelerated adjustment and prior information optimization operations is set to 0, the method of the present invention degenerates into Traditional algebraic reconstruction method. After adding the accelerated adjustment and prior information optimization operations of step (2.2.2), the algorithm speed is improved, and high-quality images are reconstructed under the condition of sparse or incomplete projection data.

First, in step (2.2.1), after performing the projection operation [x', d']=proj(x, d, ρ, a _i , b _i ) on the current vector x to the hyperplane i, there are:

||x ^* -x|| ² -||x ^* -x'|| ² =dd' (12)

The proof is as follows:

When the relaxation parameter ρ=1 in the formula (2), the vector x′ (represented as x _p ) obtained by the formula (2) is the projection point of the current vector x on the hyperplane i:

That is, x _p satisfies:

<xx _p ,x ^* -x _p >=0 (14)

Corresponding to:

||x ^* -x|| ² = ||x ^* -x _p || ² +||xx _p || ² (15)

||x ^* -x'|| ² = ||x ^* -x _p || ² +||x'-x _p || ² (16)

From formula (2) and formula (13), we can get:

Combining equations (15), (16), (17), (18), we can get:

From formula (19) and formula (3), we can get:

||x ^* -x|| ² -||x ^* -x'|| ² =d-d'

The proof is complete.

In the iterative process of the algorithm, the vector x corresponds to the auxiliary scalar d, and the auxiliary vector z corresponds to the auxiliary scalar 0. Combining the above properties, it can be seen that the following conditions are met before the accelerated adjustment and prior information optimization operations are performed in step (2.2.2):

||x ^* -z|| ² -||x ^* -x|| ² = 0-d (20)

In the following, if γ=1, the acceleration adjustment and prior information optimization operation [x', z', d']=lineAcc(x, z, d, γ) is performed in step (2.2.2) in formula (6) After that, the returned x″ is the closest point to the real solution x ^* of the equation system on the connection line between the vector x and the corresponding auxiliary vector z, that is, it satisfies

<u,x ^* -x">=0 (21)

u, x″ in the above formula are defined by formulas (4) and (6) respectively (let γ=1), that is, when γ=1, if the x″ returned by formula (6) for updating x is the vector x and Corresponding to the point on the connection line of the auxiliary vector z that is closest to the real solution x ^* of the equation system, then β in equation (6) is defined by equation (5).

The proof is as follows:

Combining formulas (21), (4), (6), we can get:

||x ^* -x|| ² =||x ^* -x″|| ² +β ² ||u|| ² (22)

||x ^* -z|| ² =||x ^* -x″|| ² +(1+β) ² ||u|| ² (23)

From equations (22), (23), we can know that:

||x ^* -z|| ² -||x ^* -x|| ² =(2β+1)||u|| ² (24)

Combining formulas (20) and (24), it can be known that:

It is the calculation expression of formula (5).

The proof is complete.

Considering the reason of noise, relaxation parameters ρ and γ are introduced into equations (2) and (6). When 0<γ≤1, the vector x″ returned by the acceleration adjustment and prior information optimization operation [x', z', d']=lineAcc(x, z, d, γ) formula (6) is smaller than before the update The vector x of is closer to the true solution x ^* of the system of equations.

The beneficial effects of adopting this technical solution:

1. Compared with the prior art (that is, the method disclosed by the inventor's previous patent application (application number 201510988996.6)), the method of the present invention directly performs multiple acceleration adjustments and priori information optimization operations within one loop iteration, due to the iterative process. It does not depend on the previous iteration variables, and it is easier to directly combine the prior information to achieve fast image reconstruction in the case of incomplete projection.

2. Compared with the prior art (that is, the method disclosed by the inventor's previous patent application (application number 201510988996.6)), the method of the present invention appropriately introduces relaxation parameters to the noise environment and derives the corresponding formulas for auxiliary scalar update and accelerated adjustment operations, Expanded the scope of application of the method.

3. Improve the convergence speed of the image reconstruction algorithm: compared with the traditional algebraic reconstruction method, the method of the present invention only adds an auxiliary vector z, and the additional storage space added can be completely ignored. In practical applications, the number T of the internal acceleration adjustment and prior information optimization operations in each iteration can be set very small (for example, less than 10), so that the additional calculation amount can also be ignored, and the simulation results show that the convergence speed of the algorithm is not can be significantly improved.

4. The application of the imaging system can be effectively expanded: the method of the present invention can be applied to different linear imaging systems, flexibly utilize the corresponding prior information, and realize fast image reconstruction under the condition of incomplete data.

Description of drawings

The present invention will be further described in detail below in conjunction with the accompanying drawings.

1 is a schematic flowchart of a fast image reconstruction method for a linear imaging system according to the present invention;

2 is a schematic flowchart of performing an update operation on a current vector x for hyperplane i=s _l in the method of the present invention;

3 is a comparison diagram of a reconstructed image and an original image according to an embodiment of the method of the present invention;

FIG. 4 is an embodiment of the method of the present invention and an algorithm square vector error convergence curve diagram combining ART and TV norm optimization.

Detailed ways

In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described below with reference to the accompanying drawings and specific embodiments.

In this embodiment, as shown in FIG. 1 and FIG. 2, a fast image reconstruction method for CT imaging system is firstly established to establish a discrete mathematical model and a priori information mathematical model of the CT imaging system, and transform the CT imaging problem into Solve large-scale linear equations and prior information optimization problems, then combine accelerated adjustment operators and prior information optimization operators to iteratively solve large-scale linear equations, optimize prior information, and reconstruct under the condition of sparse or incomplete projection data For high-quality images, the large-scale linear equations to be solved are expressed as Ax=b, where x is the vector representation of the image to be solved, which is an N×1 vector, obtained by rearranging the two-dimensional image in one dimension, and A is M ×N system matrix, the vector a _i represents the transpose of the ith row vector of matrix A, b is the M×1 measurement data vector, and b _i represents the ith element of vector b; the solution process includes the following steps:

Step 1: Initialization: Set the initial solution vector x ⁰ , relaxation parameters ρ, γ, and the number T of acceleration adjustment and prior information optimization operations performed in each loop iteration, and set the number of iterations k=0, the current vector x =x ⁰ , auxiliary scalar d=0, auxiliary vector z=x ⁰ ;

In the first step, the value range of the relaxation parameters ρ and γ is (0, 1], and the value range of the number T of acceleration adjustment and prior information optimization operations is 1≤T≤M.

Step 2: Iterative process: Repeat the following steps for the k-th iteration until the convergence conditions are met:

(2.1) Set the hyperplane projection order, and use s _l to denote the serial number of the l-th projection hyperplane, where l=1, 2,...M, M denotes the number of hyperplanes, and let t=1, t is used to represent the t-th time to perform accelerated adjustment and prior information optimization operations within the iteration;

In the step (2.1), a fixed projection sequence or a random projection sequence is used to set the hyperplane projection sequence.

(2.2) For each hyperplane i=s _l , l=1, 2,...M, perform the following operations to update the current vector x in sequence:

(2.2.1) Perform the projection operation on the current vector x to the hyperplane i, and update the current vector x and the corresponding auxiliary scalar d:

[x,d]←proj(x,d,ρ,a _i ,b _i )

Among them, [x', d']=proj(x, d, ρ, a _i , b _i ), represents the update operation of projection and auxiliary scalar d, and "←" is the assignment symbol;

The projection operation [x', d']=proj(x, d, ρ, a _i , b _i ) in the step (2.2.1) is defined as follows:

[x',d']=proj(x,d,ρ,a _i ,b _i )

λ=b _i -<a _i ,x> (1)

When ρ=1, the output x' of the projection operation [x', d']=proj(x, d, λ, a _i , b _i ) is the projection point of the current vector x to the hyperplane i. After performing the projection operation, satisfy: ||x ^* -x|| ² -||x ^* -x'|| ² =d-d', where x ^* is the true solution of the system of equations.

则对当前向量x进一步执行加速调整与先验信息优化操作，更新当前向量x，及对应的辅助向量z与辅助标量d：(2.2.2) If

Then further perform acceleration adjustment and prior information optimization operations on the current vector x, update the current vector x, and the corresponding auxiliary vector z and auxiliary scalar d:

[x,z,d]←lineAcc(x,z,d,γ)

t←t+1;

表示向上取整；[x',z',d']＝lineAcc(x,z,d,γ)表示加速调整与先验信息优化操作；In the above formula,

Represents rounding up; [x', z', d']=lineAcc(x, z, d, γ) represents accelerated adjustment and prior information optimization operations;

The acceleration adjustment and prior information optimization operation [x', z', d']=lineAcc(x, z, d, γ) in the 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)

Accelerated adjustment and prior information optimization operation [x', z', d']=lineAcc(x, z, d, γ) Formulas (4), (5), (6) perform the accelerated adjustment operation, when || x ^* -z|| ² -||x ^* -x|| ² =0-d (this formula is always satisfied before the acceleration adjustment and prior information optimization operation is performed in the iterative process of the algorithm) and when γ=1, formula (6 ) obtained x″ is the closest point to the real solution x ^* of the system of equations on the connection line between the vector x and the corresponding auxiliary vector z; when 0<γ≤1, x″ is closer than the vector x before the update to the real solution x of the system of equations ^* is closer. Formula (7) of [x', z', d']=lineAcc(x, z, d, γ) performs a priori information optimization operation, and x' returned by x'=Sup(x") is generally better than x" and The prior information of the image is closer. The prior information of an image has various representations, and is often described by some mathematical transformation, such as the 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, that is, the returned x' has a smaller value than x" Image TV norm. Formulas (8) and (9) of [x', z', d']=lineAcc(x, z, d, γ) update the auxiliary vector z and the auxiliary scalar d, respectively. After performing the accelerated adjustment and prior information optimization operation [x', z', d'] = lineAcc(x, z, d, γ), in general, the updated vector x' has more Small measurement error, which is closer to the prior information of the image (eg, the TV norm is smaller).

The step (2.1) adopts a fixed projection sequence to set the hyperplane projection sequence, and the value of the serial number i=s _l of the l (l=1, 2,...M)th projection hyperplane in each iteration is fixed, Then the acceleration adjustment and prior information optimization operations in the following step (2.2.2) in each iteration are aimed at the same hyperplane number; otherwise, if the random projection order is used to set the hyperplane projection order, the lth in each iteration The sequence number _i =sl of the projection hyperplane is random, and the acceleration adjustment and prior information optimization operations in the following step (2.2.2) in each iteration are aimed at different hyperplane numbers.

(2.3) k←k+1, that is, the number of iterations is increased by 1, and the next iteration is ready.

The present invention will be further described below with reference to a specific example.

This embodiment corresponds to the situation where the projection data is incomplete. The CT image to be reconstructed is a 128×128 two-dimensional image, that is, the original image in FIG. 3 , which corresponds to a 16384×1 vector x ^* , and A is a dimension of 7200× The system matrix of 16384, that is, M=7200, N=16384, b=Ax ^* +e is a 16384×1 measurement data vector (matrix A contains a small number of all-zero rows, which will be eliminated during processing). e is a zero-mean white Gaussian noise vector with an amplitude intensity of 1% of the measured data b.

In the case of incomplete projection, in order to reconstruct high-quality images, it is also necessary to introduce some prior knowledge for regularization. In this example, total variation (TV) regularization technique is used, and the CT image reconstruction problem after regularization is transformed into the minimization problem of the following objective function:

In this embodiment, η=0.004, the parameters of the method of the present invention in the initialization stage are set as: ρ=1, γ=1, the initial vector x ⁰ is set as an all-zero vector, and the acceleration adjustment and a priori performed in each loop iteration The number of information optimization operations is T=4, and a random projection sequence is used.

In each iteration process, a group of random permutations and combinations of 1, 2, ... M are generated as the projection order, and then the projection and acceleration adjustment and prior information optimization operations are performed on each hyperplane according to the randomly generated projection order.

)以外的投影超平面，只执行投影操作(2.2.1)，不执行加速调整与先验信息优化操作(2.2.2)；对于第1800、3600、5400、7200投影超平面，在执行投影操作(2.2.1)后还执行加速调整与先验信息优化操作(2.2.2)。由于采用的是随机投影顺序，每次迭代第1800、3600、5400、7200投影超平面的序号一般不同。比如，在第一次迭代第1800个投影超平面的序号为353，在第二次迭代第1800个投影超平面的序号为3478，这表示在第一次迭代时对序号为353的超平面进行加速调整与先验信息优化操作，而在第二次迭代时对序号为3478的超平面进行加速调整与先验信息优化操作。For the 1800th, 3600th, 5400th, 7200th (ie

), only the projection operation (2.2.1) is performed, and the acceleration adjustment and prior information optimization operations (2.2.2) are not performed; for the 1800th, 3600, 5400, and 7200 projection hyperplanes, the projection operation is performed After (2.2.1), the accelerated adjustment and prior information optimization operations (2.2.2) are also performed. Due to the random projection order, the serial numbers of the 1800th, 3600th, 5400th, and 7200th projection hyperplanes are generally different in each iteration. For example, the sequence number of the 1800th projection hyperplane in the first iteration is 353, and the sequence number of the 1800th projection hyperplane in the second iteration is 3478, which means that in the first iteration, the hyperplane with sequence number 353 is processed. Accelerate adjustment and prior information optimization operations, and perform accelerated adjustment and prior information optimization operations on the hyperplane with serial number 3478 in the second iteration.

Since the method of the present invention does not depend on the previous iteration variable in the iterative process, it is easy to directly combine the prior information to realize the fast image reconstruction under the condition of incomplete projection. The comparison diagram of the image reconstructed by the method of the present invention and the original image is shown in Figure 3, and the square vector error convergence curve diagram of the method of the present invention and the algorithm combining ART and TV norm optimization is shown in Figure 4. It can be seen from the figure The method of the present invention significantly improves the convergence speed.

The foregoing has shown and described the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above-mentioned embodiments, and the descriptions in the above-mentioned embodiments and the description are only to illustrate the principle of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have Various changes and modifications fall within the scope of the claimed invention. The claimed scope of the present invention is defined by the appended claims and their equivalents.

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 dimension_iRepresenting the transpose of the ith row vector of matrix A, b being the M × 1 measurement data vector, using b_iThe ith element representing vector b; the solving process comprises the following steps:

the method comprises the following steps: initialization: setting an initial solution vector x⁰Relaxation 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 x⁰Auxiliary scalar d is 0 and auxiliary vector z is x⁰；

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 s_lA 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 ═ s_lM, 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,ρ,a_i,b_i)

wherein, [ x ', d']＝proj(x,d,ρ,a_i,b_i) 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 iteration_lIf 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 s_lIs 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,ρ,a_i,b_i) The definition is as follows:

[x',d']＝proj(x,d,ρ,a_i,b_i)

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.

Images