CN105590332A - Rapid algebraic reconstruction technique applied to computed tomography imaging - Google Patents

Abstract

The invention provides a rapid algebraic reconstruction technique (ART) applied to computed tomography (CT) imaging, belonging to the technical field of biomedical imaging, nondestructive test, etc. Based on traditional ART, the rapid ART of the invention randomly selects partial hyperplane, and performs further acceleration adjustment operation on the projection solution vector obtained through the traditional ART on the partial hyperplane, thereby allowing the solution vector after adjustment and update to be optimal on the connection line between an original solution vector and the hyperplane solution vector of last iteration. According to the rapid ART, a different hyperplane projection sequence can be employed for each iteration. A random projection sequence of each iteration allows higher efficiency under various conditions. Compared with the traditional ART, the rapid ART substantially increases an algorithm convergence speed, and can be applied to rapid CT image reconstruction under incomplete projection conditions.

Description

A kind of quick algebraic reconstruction method that is applied to CT imaging

Technical field

The present invention relates to a kind of quick algebraic reconstruction method of the CT of being applied to imaging, belong to biomedical imaging,The technical fields such as Non-Destructive Testing.

Background technology

Computed tomography (CT) imaging is as a kind of important Dynamic Non-Destruction Measurement, in medical domain and workIndustry field is all widely used. Conventional CT image reconstruction algorithm mainly contains two large classes at present: a class is to filterThe analytic method that ripple backprojection algorithm (FBP) is representative; Another kind of taking algebraic reconstruction technique (ART) as representativeIterative algorithm. The advantage of analytic method is that amount of calculation is little, and reconstruction speed is fast, can obtain good in the time that data for projection is sufficientGood reconstructed image quality, therefore main flow business CT still adopts and resolves class methods reconstruction image now.

In practical medical application and commercial Application, the factors such as examined environment, cost, time, personnel healthLimit, often run into the situation of incomplete projection data. In incomplete projection situation, existing filtering is counter throwsShadow class algorithm is difficult to reconstruct complete CT image. Different from analytic method, algebraic reconstruction algorithm is by image reconstructionProblem is converted into extensive Solving Linear problem, by conjunction with some regularization prioris, in projectionWhen data transformation, still can rebuild high quality graphic. In addition, algebraic reconstruction algorithm is not subject to reason problem model of living inRestriction and can introduce the prior information of tested object, be easily extended to the other field of CT application, exhibitionShow very large development space and application potential.

The shortcoming of algebraic reconstruction technique is that amount of calculation is large, and reconstruction speed is slow, and this has limited it within a very long timeApplication. But along with the development of parallel computation theory and computer hardware technology, this type of technical bottleneck has obtained veryThe alleviation of large degree, the advantage of algebraic reconstruction technique highlights more, and algebraic reconstruction technique causes numerous again in recent yearsScholar's attention. The convergence rate and the precision that how in incomplete projection situation, to improve algebraic reconstruction technique, becomeThe algebraic reconstruction algorithm key that can be widely used in actual imaging system.

Summary of the invention

The invention provides a kind of quick algebraic reconstruction method of the CT of can be applicable to imaging. With respect to traditional algebraicallyMethod for reconstructing (ART), this algorithm the convergence speed is fast, can be for the rapid CT figure in incomplete projection situationPicture is rebuild.

Theory analysis of the present invention is:

Image reconstruction problem is converted into following extensive Solving Linear problem by algebraic reconstruction technique (ART):

Ax＝b(1)

Wherein, x is image to be asked, and is N × 1 vector, by two dimension or 3-D view rearranged by one dimensionArrive. Each element x of x_iA pixel of corresponding original image. A is M × N sytem matrix, each element a_ijJ the contribution of pixel to i article of ray described. B is M × 1 measurement data vector, its i element b_iThe decay of corresponding i article of ray.

ART method is the earliest proposed by Kaczmarz, also claims Kaczmarz method (being called for short KA). It is basicThought is: each equation is regarded as to a hyperplane of N dimension space, from initial solution vector x⁰Set out, exist successivelyProjection on each hyperplane, progressively approaches the true solution of equation group, and available following mathematical formulae represents:

x^k+1＝P(b,x^k)＝P_Mο…οP₂οP₁(b,x^k)(2)

x^k,i＝P_i(b,x^k,j)＝x^k,j+λ^i,ja_i(3)

${\mathrm{\&lambda;}}^{i,j}=\frac{b_i-<a_i,x^{k,j}>}{\left|\right|a_i|{|}^2}=\frac{b_i-a_i^T{x}^{k,j}}{\left|\right|a_i|{|}^2}---\left(4\right)$

x^k，0＝x^k,x^k+1＝x^k,M(5)

Wherein, a_iFor the i row vector of matrix A,<>with || || inner product and Euler's norm represented respectively, a_i ^TIt is right to representa_iTransposition. P_i(b,x^k,j) represent vector x^k,jProject to the operation of i hyperplane, x^k,jRepresent KA kWhen inferior iteration, project to the vector of j hyperplane. Each iteration, to vector x^kExecutable operations P (b, x^k),On each hyperplane, carry out taking turns successively projection, obtain follow-on solution vector x^k+1＝P(b,x^k)＝x^k,M。

In order to improve the efficiency of Kaczmarz method, it is individual super flat that the present invention projects to i in each iteration to KAThe projection vector x of face^k,iAlong itself and this hyperplane projection vector of former generation x^k-1,iLine direction is further adjustedOptimize, make the new vector x after adjusting^cAt x^k,i-1With x^k,iOptimum (to true solution vector x in line direction^*Squared-distance or error minimum). True solution vector x^*For the intersection point of all hyperplane, adjustment obtains after optimizingVector x^cShould meet: (x*-x^c)⊥(x^k,i-x^k-1,i). Geometric representation is as shown in accompanying drawing Fig. 1.

The iteration mathematical formulae that the present invention puies forward quick algebraic reconstruction method is described below:

$x^c={\tilde{P}}_i^k(b,x^{k,j})=P_i(b,x^{k,j})+{\mathrm{\&beta;x}}^e---\left(7\right)$

x^e＝P_i(b,x^k,j)-x^k-1,i＝x^k,i-x^k-1,i(8)

The key of this algorithm is how to try to achieve optimum β. Utilize x^cAbove-mentioned vertical relation, have:

d^k,i＝||x*-x^k,i||²＝||x*-x^c||²+β²||x^e||²(9)

d^k-1,i＝||x*-x^k-1,i||²＝||x*-x^c||²+(β+1)²||x^e||²(10)

△d＝d^k-1,i-d^k,i＝(2β+1)||x^e||²(11)

$\mathrm{\&beta;}=\frac{\mathrm{\&Delta;}d-\left|\right|x_e|{|}^2}{2\left|\right|x_e|{|}^2}---\left(12\right)$

Wherein, d^k,iRepresent vector x^k,iTo true solution vector x^*Square error (or distance). According to formula (3), (4),In KA by vector x^k,jWhile projecting to i plane, the x obtaining^k,iMeet: (x^k,i-x^k,j)⊥(x^k,i-x*)，Therefore have:

d^k,i＝d^k,j-||x^k,i-x^k,j||²＝d^k,j-||λ^i,ja_i||²(13)

Like this by an initial value d⁰Start (d^k,0＝d^k), according to formula (13) iteration successively, can obtain adjacentTwice iteration d^k,iWith respect to d^k-1,iReduce value:

△d＝d^k-1,i-d^k,i(14)

In formula (14), d⁰Can be eliminated. Although so d⁰Value the unknown, but can be by d in the time that programming realizes⁰IfFor any one constant C, the result of formula (14) can not change, and does not also affect the iterative computation of algorithm. WillThe △ d substitution formula (12) that formula (14) calculates can be obtained β, then the β substitution formula (7) of obtaining is obtainedx^c. Then upgrade x according to formula (15), (16)^k,iWith d^k,i：

x^k,i＝x^c(15)

d^k,i＝d^k,i-β²||x^e||²(16)

It is pointed out that the d that generally programming is set⁰≠||x*-x⁰||², algorithm iteration obtains like thisd^k,iWith || x*-x^k,i||²Differ a constant, but △ d=d^k-1,i-d^k,i＝||x*-x^k-1,i||²-||x*-x^k,i||²Set up.

In the quick algebraic reconstruction technique that the present invention carries, each iteration can adopt different hyperplane projection sequence.Emulation and initial analysis show, adopt random projection sequence in each iteration, in many cases can be obviousImprove the efficiency of algorithm.

The algorithm that the present invention carries, each hyperplane i need to record the projection vector x of a front iteration^k-1,i. RightIn CT imaging applications, the size of image (is vector x^k-1,iDimension) and the number of hyperplane (measure numberAccording to number) very large. In this case, record the front once iterative vectorized x of each hyperplane^k-1,iNeed very largeMemory space. In order to reduce memory requirement, can promote this algorithm, only select a small amount of hyperplane to aboveProjection vector x^k,iCarry out further adjusting and optimizing. Can accelerate degree and required extra storage sky at algorithm like thisBetween carry out balance, accelerate convergence of algorithm speed increasing under the condition of a small amount of memory space. Fortunately, rightIn algorithm of the present invention, this mode is easy to realize, and only needs to select in advance one and accelerates to adjust the super of operationPlane subset Ω, iteration more new formula is constant.

The technical scheme of the quick algebraic reconstruction method that the present invention carries comprises the following steps:

Step 1: initialize. Set k=0, set initial solution x^0,0＝x⁰，d^0,0＝d⁰, arrange and accelerate to adjustThe hyperplane subset Ω of operation.

Step 2: iterative process. Repeat following steps until algorithm meets the condition of convergence:

1. the k time iteration, first set hyperplane projection sequence, represent that with idx (l) l projection is super flatThe index subscript of face.

2. the k time iteration, to each hyperplane i=idx (l), l=1,2 ... M, carries out successively as finishes drillingDo:

1) by current vector x^k,jProject to hyperplane i, obtain vector x^k,i. Concrete operations are as follows:

a)j＝idx(l-1)

b) $\mathrm{\&lambda;}=\frac{b_i-<a_i,x^{k,j}>}{\left|\right|a_i|{|}^2}$

c)x^k,i＝x^k,j+λa_i

d)d^k,i＝d^k,j-||λa_i||²

2) if k > 0 and i ∈ Ω, to the vector x on hyperplane i^k,iCarry out the following operation of adjusting:

a)△d＝d^k-1,i-d^k,i

b)x^e＝x^k,i-x^k-1,i

c) $\mathrm{\&beta;}=\frac{\mathrm{\&Delta;}d-\left|\right|x_e|{|}^2}{2\left|\right|x_e|{|}^2}$

d)x^k,i＝x^k,i+βx^e

e)d^k,i＝d^k,i-β²||x^e||²

3. upgrade vector, iterations adds 1. Concrete operations are as follows:

x^k+1＝x^k,M，d^k+1＝d^k,M

x^k+1，0＝x^k+1，d^k+1，0＝d^k+1

k＝k+1

The present invention can select less hyperplane to carry out and accelerate to adjust operation its projection vector, and at every turnIteration adopts accidental projection order, can reduce like this amount of extra memory, shortens convergence time.

Brief description of the drawings

Fig. 1 projection vector is adjusted geometric representation.

Square vector distance (or error) of table 1 iterative process

Detailed description of the invention

Below in conjunction with the drawings and specific embodiments, the present invention is further detailed.

Suppose that it is Ax=b that CT imaging problem transforms the extensive system of linear equations obtaining. Wherein, A is M × NSytem matrix, the dimension that M is measurement data, the dimension that N=n × n is image space. In incomplete projection situationUnder, M < N. In order to rebuild high quality graphic, also need to introduce some prioris and carry out regularization. One hasEffect and widely used Regularization Technique are full variation (Totalvariance, TV) Regularization Techniques.

The invention provides the Regularization Technique of a kind of similar TV, after regularization CT image reconstruction problem be converted into asThe minimization problem of lower object function:

f(x)＝||Ax-b||²+γ²||R₁x||²+γ²||R₂x||²(17)

Wherein, γ is regularization coefficient. R₁、R₂For N × N matrix. Suppose that former n × n two dimensional image is by rearrangementBecome N × 1 vector x, make h_r(i) represent i adjacent pixel subscript in pixel the right of original image, v_b(i) represent theI the adjacent pixel subscript in pixel below. Matrix R₁Element definition as follows:

$r_{i,j}=\{\begin{array}{cc}1,& ifj=h_r\left(i\right)\\ -1,& ifj=j\\ 0,& else\end{array},i\&le;N-n---\left(18\right)$

Matrix R₂Element definition as follows:

${r^{\&prime;}}_{i,j}=\left\{\begin{array}{cc}1,& ifj=v_b\left(i\right)\\ -1,& ifj=i\\ 0,& else\end{array}\right.,\mathrm{mod}(i,n)\&NotEqual;0---\left(19\right)$

In above formula, mod (i, n) represents the remainder of i divided by n.

Object function (17) is carried out to differentiate to vector x:

▽f(x)＝2(A^TA+γ²R₁ ^TR₁+γ²R₂ ^TR₂)x-2A^Tb(20)

The vector x that gradient ▽ f (x) is 0 is the optimal solution that correspondence minimizes object function (17). Minimize object function(17), be equivalent to and solve following system of linear equations:

Cx＝q(21)

Wherein,

C＝(A^TA+γ²R₁ ^TR₁+γ²R₂ ^TR₂)(22)

q＝A^Tb(23)

The detailed description of the invention of algorithm of the present invention is described with a simple example below. Suppose to obtain after regularizationEquation group (21) in C be following 3 × 3 matrixes:

$C=\left[\begin{array}{ccc}1& 2& 5\\ 2& 7& 9\\ 3& 8& 6\end{array}\right]$

Optimal vector x*=[123]^T，q＝Cx*＝[204337]^T。

The quick algebraic reconstruction method that adopts the present invention to propose solves above-mentioned equation group, comprises following step:Step 1: initialize. Set k=0, initial solution x^0,0＝x⁰＝0，d^0,0＝d⁰=0, accelerate to adjustThe hyperplane subset Ω={ 2} of operation.

Step 2: iterative process. Repeat following steps until algorithm meets the condition of convergence:

1. each iteration, first sets hyperplane projection sequence. Hyperplane projection sequence in this example is fixedly adoptedWith pressing i=1,2,3 ascending order arrangement modes.

2. the 1st iteration (k=0). To each hyperplane i=1,2,3, carry out successively following operation:

1) by x^0,0=0 projects to hyperplane i=1, the projection vector x obtaining^0,1And d^0,1Be respectively:

x^0,1＝[0.6671.33333.3333]^T

d^0,1＝-13.3333

Iteration for the first time, does not accelerate to adjust operation to projection vector.

2) by x^0,1Project to hyperplane i=2, the projection vector x obtaining^0,2And d^0,2Be respectively:

x^0,2＝[0.70151.45523.4900]^T

d^0,2＝-13.374

3) by x^0,2Project to hyperplane i=3, the projection vector x obtaining^0,3And d^0,3Be respectively:

x^0,3＝[0.76521.62503.6174]^T

d^0,3＝-13.4231

4) upgrade vector, iterations adds 1. Concrete operations are as follows:

x¹＝x^1,0＝x^0,3＝[0.76521.62503.6174]^T

d¹＝d^1,0＝d^0,3＝-13.4231

k＝k+1

3. the 2nd iteration (k=1). To each hyperplane i=1,2,3, carry out successively following operation:

1) by x^1,0Project to hyperplane i=1, the projection vector x obtaining^1,1And d^1,1Be respectively:

x^1,1＝[0.69511.48493.2670]^T

d^1,1＝-13.5704

Because i=1 does not belong to, hyperplane subset Ω={ 2}, therefore to projection vector x^1,1AddThe whole operation of velocity modulation.

2) by x^1，1Project to hyperplane i=2, the projection vector x obtaining^1,2And d^1,2Be respectively:

x^1,2＝[0.72211.57963.3888]^T

d^1,2＝-13.5949

Because i=2 belongs to, hyperplane subset Ω={ 2}, therefore then will be to projection vector x^1,2Carry outAccelerate to adjust operation, adjust the vector x after upgrading^1,2And d^1,2Be respectively:

x^1,2＝[0.79912.04273.0114]^T

d^1,2＝-13.9577

3) by x^1,2Project to hyperplane i=3, the projection vector x obtaining^1,3And d^1，3Be respectively:

x^1,3＝[0.80442.05683.0220]^T

d^1,3＝-13.9580

Because i=3 does not belong to hyperplane subset Ω, therefore not to projection vector x^1,3Accelerate to adjustOperation.

4) upgrade vector, iterations adds 1. Concrete operations are as follows:

x²＝x^2,0＝x^1,3＝[0.80442.05683.0220]^T

d²＝d^2,0＝d^1,3＝-13.9580

k＝k+1

4. the 3rd iteration (k=2)

……

In this example, initial solution vector x⁰=0 to true solution vector x*=[123]^TSquared-distance be:

D⁰＝||x*||²＝14

D^k,i＝d^k,i+ 14 are x^k,iTo the squared-distance (or error) of true solution vector x*. Accompanying drawing table 1 is listedFor the first time with the each hyperplane solution vector of iteration x for the second time^k,iTo the squared-distance D of true solution vector x*^k,i。As can be seen from the table, from initial vector x^0,0＝x⁰=0 starts the mistake at each hyperplane successive projectionCheng Zhong, D^k,iConstantly reducing. D^0,1=0.6667, take turns after iteration D through one^1,1＝0.4296。x^1,1Project to the projection vector x that the 2nd hyperplane obtains^1,2Corresponding D^1,2=0.4051, reduce to D^1,1's94%. Adopt algorithm of the present invention to x^1,2After accelerating to adjust, the x of renewal^1,2Corresponding D^1,2＝0.0423，Square vector error significantly reduces, and reduces to 10% before renewal. Visible, the method that the present invention carries is passableObviously improve the speed of algorithm for reconstructing.

Table 1

Claims (3)

1. be applied to a quick algebraic reconstruction method for CT imaging, i.e. the solution of iterative equation group Ax=b, wherein, x isFor image to be asked, be N × 1 vector, by two dimension or 3-D view are rearranged and obtained by one dimension; Each element x of x_iCorrespondingA pixel of original image; A is M × N sytem matrix, each element a_ijJ the contribution of pixel to i article of ray described; bFor M × 1 measurement data vector, its i element b_iThe decay of corresponding i article of ray; The method comprises:

Step 1: initialize: set k=0, set initial solution x^0,0＝x⁰，d^0,0＝d⁰, arrange and accelerate to adjust the super of operationPlane set omega;

Step 2: iterative process: repeat following steps until meet the condition of convergence:

Step 2.1, the k time iteration, is first set hyperplane projection sequence, represents the index of l super projection plane with idx (l)Subscript;

Step 2.2 is the k time iteration, to each hyperplane i=idx (l), and l=1,2 ... M, carries out following operation successively:

Step 2.2.1 is by current solution vector x^k,jProject to hyperplane i, obtain solution vector x^k,i; Concrete operations are as follows:

A) make under the index of l-1 super projection plane and be designated as j;

B) adopt formulaCalculate λ, a_iFor the i row vector of matrix A;

C) adopt formula x^k,i＝x^k,j+λa_iCalculate solution vector x^k，i；

D) calculate d^k,i＝d^k,j-||λa_i||², wherein, d^k,iRepresent solution vector x^k,iTo true solution vector x^*Square error;

If step 2.2.2 is k > 0 and i ∈ Ω, to the vector x on hyperplane i^k,iCarry out the following operation of adjusting:

A) calculate △ d=d^k-1,i-d^k,i；

B) calculate x^e＝x^k,i-x^k-1,i；

C) calculate $\mathrm{\&beta;}=\frac{\mathrm{\&Delta;}d-\left|\right|x_e|{|}^2}{2\left|\right|x_e|{|}^2};$

D) upgrade x^k,i＝x^k,i+βx^e；

E) upgrade d^k，i＝d^k,i-β²||x^e||²；

Step 2.3 is upgraded vector, and iterations adds 1;

x^k+1＝x^k,M，d^k+1＝d^k,M

x^k+1，0＝x^k+1，d^k+1，0＝d^k+1

k＝k+1。

2. a kind of quick algebraic reconstruction method that is applied to CT imaging as claimed in claim 1, is characterized in that in step 2.1Adopt the projection sequence of first determining all hyperplane, then according to this sequence, each hyperplane is carried out to projection, wherein determine all superThe method of plane projection order is permanent order or random sequence.

3. a kind of quick algebraic reconstruction method that is applied to CT imaging as claimed in claim 1, is characterized in that establishing in step 1The hyperplane set omega of putting is any nonvoid subset of hyperplane complete or collected works.