CN105701847A - Algebraic reconstruction method of improved weight coefficient matrix - Google Patents

Abstract

The invention relates to an algebraic reconstruction method of an improved weight coefficient matrix and belongs to the CT projection data reconstruction method technology field. The method mainly aims at solving problems of logical loopholes in a Ray-Box Intersection weight coefficient matrix algorithm and a fuzzy SART+TVM reconstruction image detail based on a bilinear interpolation. Through improving the Ray-Box Intersection weight coefficient matrix algorithm, reconstruction time is shortened, the image detail is maintained and a signal to noise ratio is increased. Simultaneously, a TVM algorithm is used as a bound term, an iteration reconstruction convergence speed is accelerated, artifacts and noises of the reconstruction image are reduced and reconstruction image quality is effectively increased.

Description

A kind of arithmetic reconstruction method improving weight coefficient matrix

Technical field

The invention belongs to CT data for projection method for reconstructing technical field, relate to a kind of associating arithmetic reconstruction method based on the Ray-BoxIntersection weight coefficient matrix improved。

Background technology

Computerized tomography (ComputedTomography, CT) technology is a kind of important Dynamic Non-Destruction Measurement, its core is Image Reconstruction Technology, namely utilizes the data for projection under multiple sampling visual angle to reconstruct two dimension or the 3-D view of object inside and outside portion architectural feature。CT algorithm for reconstructing mainly includes analytic reconstruction algorithm and iterative reconstruction algorithm。Analytic reconstruction needs the data for projection that uniform sampling is enough within the scope of 180 ° or 360 °, but in industrial CT system practical application, often occur that workpiece is oversize makes sampling angle limited owing to scanning, or X ray cannot penetrate high-density region and cause that sample count is not enough, or reduce sampling density to improve the situations such as scanning speed, these situations belong to incomplete data for projection Problems of Reconstruction, it is impossible to meet the condition of analytic reconstruction, need to adopt iterative reconstruction algorithm。

In recent years, solve incomplete data for projection Problems of Reconstruction, it is generally adopted the iterative reconstruction algorithm based on image total variation, it is based especially on associating algebraic reconstruction algorithm (the SimultaneousAlgebraicReconstructionTechniqueBasedTotalVa riationMinimization of image total variation minimization, SART+TVM), it is openness that this algorithm utilizes that CT image gradient converts, and combine with SART algorithm, in each iterative process, preliminary reconstruction image is obtained first with SART algorithm, then TVM algorithm is utilized to adjust image total variation along image gradient direction。SART+TVM algorithm reconstructed image quality has multiple influence factor, including: SART iterations, SART relaxation factor, weight coefficient matrix, TVM iterations, TVM regulatory factor etc., and weight coefficient matrix is one of key factor determining SART+TVM algorithm reconstructed image quality。

Bilinear interpolation algorithm is the presently most used method solving weight coefficient matrix, it utilizes neighbor pixel to carry out interpolation by proportional parts and solves the pixel value of current point, the more complete structural information of ratio can be obtained, but image border easily obscures, by contrast, Ray-BoxIntersection algorithm in principle closer to real projection process, but because it exists logic leak, such as: ray is non-intersect with certain block of pixels, but it is judged as intersecting, or weighted value is more than the actual (tube) length angle value of ray traverse block of pixels, causes that reconstructed image quality is poor, artifact is many。Accordingly, it would be desirable to Ray-BoxIntersection algorithm is improved, making up the logic leak himself existed, improve reconstructed image quality, especially rebuilding image detail to highlight。

Summary of the invention

In view of this, it is an object of the invention to provide a kind of arithmetic reconstruction method improving weight coefficient matrix, the logic leak that the method exists for Ray-BoxIntersection weight coefficient matrix algorithm and the SART+TVM based on bilinear interpolation rebuild the problem that image detail is fuzzy, Ray-BoxIntersection weight coefficient matrix algorithm is improved, thus shortening reconstruction time, keep image detail, improve signal to noise ratio；Utilize TVM algorithm as bound term simultaneously, accelerate iterative approximation convergence rate, decrease the artifact and noise of rebuilding image, be effectively increased reconstructed image quality。

For reaching above-mentioned purpose, the present invention provides following technical scheme:

A kind of arithmetic reconstruction method improving weight coefficient matrix, comprises the following steps:

S1: obtain data for projection p_i, initialization CT scan parameter, i=0,1,2 ..., N, N is projection view angles sum；

S2: by image x to be reconstructed_jCompose initial value,J=0,1,2 ..., M-1, M is total number of image pixels, and k is SART iterations；

S3: according to ray-driven mode, calculate the weight coefficient matrix A={a under this projecting direction_ij, a_ijIt it is i-th ray weighted value to jth block of pixels；

S4: orthographic projection, obtains the simulated projections value of i-th ray

S5: according to the actual measurement projection value p of ray_i, simulated projections valueWith weight coefficient matrix A, obtain the correction value of i-th ray ${\mathrm{\Δ}}_i=(p_i-{\tilde{P}}_i)/\underset{m=1}{\overset{M}{\Σ}}{a}_{im}{x}_m^{\left(k\right)};$

S6:i=i+1, repeats step S3-S5, until during the correction value of all rays completed under this projecting direction, carrying out back projection's more new images according to following SART iterative formula and obtain

$x_j^{(k+1)}=x_j^k+{\mathrm{\λ}}^k\frac{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}\[a_{ij}{\mathrm{\Δ}}_i\]}{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}{a}_{ij}}$

In formulaWithThe respectively pixel value of jth block of pixels, I in kth+1 and kth subiterations process_θFor the set of all rays under angle, θ, λ^kFor relaxation factor in kth subiterations process；

S7: repeat step S2-S6, until completing the correction of all projection angles；

S8: after image carrying out gradient descent method and adjusts image total variation, image total variation formula is:

$\frac{u_{n+1}-u_n}{\mathrm{\Δ}t}=\▿\·\frac{\▿u_n}{\left|\right|\▿u_n\left|\right|}-\mathrm{\α}(u_n-x_j^{(k+1)})$

1≤n < N_TVM,N_TVMFor the iterations of TVM, u_nFor total variation minimization image, u_nRepresentative image u_nGradient, | | | | represent L1 norm operator, α be TVM reconcile the factor；

S9: orderJudge whether to reach the condition of convergenceε is the minimum more than 0, is, turns to step S11, is not, turns to step S10；

S10: reset iterative imageRepeat step S3-S9；

S11: exit circulation, obtains rebuilding image。

Further, described step S3 specifically includes:

S31: by beam angularly discretization, utilizes ray and detector relation one to one, can obtain the positional information of every ray；

S32: utilize the algebraic relation of ray and image-region, obtains ray starting point (x in image-region₁,y₁) and terminal (x₂,y₂)；

S33: ask for the lower left corner coordinate set { (Indx, Indy) } of ray all block of pixels of process in reconstruction regions；

S34: calculate the intersection point of ray and passed block of pixels；

S35: calculate the distance between two intersection points in same block of pixels, and be considered as this block of pixels weight for this ray, distinguishingly, when radiation levels or when being perpendicular to block of pixels, weight is 1；

S36: repeat step S32-S35, until traveling through all rays。

The beneficial effects of the present invention is: method provided by the invention shortens reconstruction time, maintains image detail, improves signal to noise ratio；Utilize TVM algorithm as bound term simultaneously, accelerate iterative approximation convergence rate, decrease the artifact and noise of rebuilding image, be effectively increased reconstructed image quality。

Accompanying drawing explanation

In order to make the purpose of the present invention, technical scheme and beneficial effect clearly, the present invention provides drawings described below to illustrate:

Fig. 1 is SART+TVM algorithm for reconstructing flow chart；

Fig. 2 is algorithm flow chart before improving；

Fig. 3 is algorithm flow chart after improving；

Fig. 4 is ray and reconstruction regions schematic diagram。

Detailed description of the invention

Below in conjunction with accompanying drawing, the preferred embodiments of the present invention are described in detail。

Fig. 1 is SART+TVM algorithm for reconstructing flow chart, Fig. 2 is algorithm flow chart before improving, Fig. 3 is algorithm flow chart after improving, Fig. 4 is ray and reconstruction regions schematic diagram, as shown in the figure, present invention is generally directed to incomplete data for projection Problems of Reconstruction, adopt SART+TVM algorithm for reconstructing, on the basis of SART iterative reconstruction algorithm, namely introduce the constraint of image total variation minimization。

This method specifically includes following steps:

S1: obtain data for projection p_i, initialization CT scan parameter, i=0,1,2 ..., N, N is projection view angles sum；

S2: by image x to be reconstructed_jCompose initial value,J=0,1,2 ..., M-1, M is total number of image pixels, and k is SART iterations；

S3: according to ray-driven mode, calculate the weight coefficient matrix A={a under this projecting direction_ij, a_ijIt it is i-th ray weighted value to jth block of pixels；

S4: orthographic projection, obtains the simulated projections value of i-th ray

S5: according to the actual measurement projection value p of ray_i, simulated projections valueWith weight coefficient matrix A, obtain the correction value of i-th ray ${\mathrm{\&Delta;}}_i=(p_i-{\tilde{P}}_i)/\underset{m=1}{\overset{M}{\&Sigma;}}{a}_{im}{x}_m^{\left(k\right)};$

S6:i=i+1, repeats step S3-S5, until during the correction value of all rays completed under this projecting direction, carrying out back projection's more new images according to following SART iterative formula and obtain

$x_j^{(k+1)}=x_j^k+{\mathrm{\&lambda;}}^k\frac{\underset{i\&Subset;I_{\mathrm{\&theta;}}}{\&Sigma;}\&lsqb;a_{ij}{\mathrm{\&Delta;}}_i\&rsqb;}{\underset{i\&Subset;I_{\mathrm{\&theta;}}}{\&Sigma;}{a}_{ij}}$

In formulaWithThe respectively pixel value of jth block of pixels, I in kth+1 and kth subiterations process_θFor the set of all rays under angle, θ, λ^kFor relaxation factor in kth subiterations process；

S7: repeat step S2-S6, until completing the correction of all projection angles；

S8: after image carrying out gradient descent method and adjusts image total variation, image total variation formula is:

$\frac{u_{n+1}-u_n}{\mathrm{\&Delta;}t}=\&dtri;\&CenterDot;\frac{\&dtri;u_n}{\left|\right|\&dtri;u_n\left|\right|}-\mathrm{\&alpha;}(u_n-x_j^{(k+1)})$

1≤n < N_TVM,N_TVMFor the iterations of TVM, u_nFor total variation minimization image, u_nRepresentative image u_nGradient, | | | | represent L1 norm operator, α be TVM reconcile the factor；

S9: orderJudge whether to reach the condition of convergenceε is the minimum more than 0, is, turns to step S11, is not, turns to step S10；

S10: reset iterative imageRepeat step S3-S9；

S11: exit circulation, obtains rebuilding image。

In previous step S3, weight coefficient matrix algorithm adopts the Ray-BoxIntersection algorithm improved, and it is embodied as step and is:

S31: by beam angularly discretization, utilizes ray and detector relation one to one, can obtain the positional information of every ray；

S32: utilize the algebraic relation of ray and image-region, obtains ray starting point (x in image-region₁,y₁) and terminal (x₂,y₂)；

S33: ask for the lower left corner coordinate set { (Indx, Indy) } of ray all block of pixels of process in reconstruction regions；

S34: calculate the intersection point of ray and passed block of pixels；

S35: calculate the distance between two intersection points in same block of pixels, and be considered as this block of pixels weight for this ray, distinguishingly, when radiation levels or when being perpendicular to block of pixels, weight is 1；

S36: repeat step S32-S35, until traveling through all rays。

What finally illustrate is, preferred embodiment above is only in order to illustrate technical scheme and unrestricted, although the present invention being described in detail by above preferred embodiment, but skilled artisan would appreciate that, in the form and details it can be made various change, without departing from claims of the present invention limited range。

Claims (2)

1. the arithmetic reconstruction method improving weight coefficient matrix, it is characterised in that: comprise the following steps:

S1: obtain data for projection p_i, initialization CT scan parameter, i=0,1,2 ..., N, N is projection view angles sum；

S2: by image x to be reconstructed_jCompose initial value,J=0,1,2 ..., M-1, M is total number of image pixels, and k is SART iterations；

S3: according to ray-driven mode, calculate the weight coefficient matrix A={a under this projecting direction_ij, a_ijIt it is i-th ray weighted value to jth block of pixels；

S4: orthographic projection, obtains the simulated projections value of i-th ray

S5: according to the actual measurement projection value p of ray_i, simulated projections valueWith weight coefficient matrix A, obtain the correction value of i-th ray ${\mathrm{\&Delta;}}_i=(p_i-\widetilde{P_i})/\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_{\mathrm{im}}{x}_m^{\left(k\right)};$

S6:i=i+1, repeats step S3-S5, until during the correction value of all rays completed under this projecting direction, carrying out back projection's more new images according to following SART iterative formula and obtain

$x_j^{(k+1)}=x_j^k+{\mathrm{\&lambda;}}^k\frac{\underset{i\&Subset;I_{\mathrm{\&theta;}}}{\&Sigma;}\&lsqb;a_{ij}{\mathrm{\&Delta;}}_i\&rsqb;}{\underset{i\&Subset;I_{\mathrm{\&theta;}}}{\&Sigma;}{a}_{ij}}$

In formulaWithThe respectively pixel value of jth block of pixels, I in kth+1 and kth subiterations process_θFor the set of all rays under angle, θ, λ^kFor relaxation factor in kth subiterations process；

S7: repeat step S2-S6, until completing the correction of all projection angles；

S8: after image carrying out gradient descent method and adjusts image total variation, image total variation formula is:

$\frac{u_{n+1}-u_n}{\mathrm{\&Delta;}t}=\&dtri;\&CenterDot;\frac{\&dtri;u_n}{\left|\right|\&dtri;u_n\left|\right|}-\mathrm{\&alpha;}(u_n-x_j^{(k+1)})$

1≤n < N_TVM,N_TVMFor the iterations of TVM, u_nFor total variation minimization image,Representative image u_nGradient, | | | | represent L1 norm operator, α be TVM reconcile the factor；

S9: orderJudge whether to reach the condition of convergenceε is the minimum more than 0, is, turns to step S11, is not, turns to step S10；

S10: reset iterative imageRepeat step S3-S9；

S11: exit circulation, obtains rebuilding image。

2. a kind of arithmetic reconstruction method improving weight coefficient matrix according to claim 1, it is characterised in that: described step S3 specifically includes:

S31: by beam angularly discretization, utilizes ray and detector relation one to one, can obtain the positional information of every ray；

S32: utilize the algebraic relation of ray and image-region, obtains ray starting point (x in image-region₁,y₁) and terminal (x₂,y₂)；

S33: ask for the lower left corner coordinate set { (Indx, Indy) } of ray all block of pixels of process in reconstruction regions；

S34: calculate the intersection point of ray and passed block of pixels；

S35: calculate the distance between two intersection points in same block of pixels, and be considered as this block of pixels weight for this ray, distinguishingly, when radiation levels or when being perpendicular to block of pixels, weight is 1；

S36: repeat step S32-S35, until traveling through all rays。