CN107292846B - Recovery method of incomplete CT projection data under circular orbit - Google Patents

Recovery method of incomplete CT projection data under circular orbit Download PDF

Info

Publication number
CN107292846B
CN107292846B CN201710501987.9A CN201710501987A CN107292846B CN 107292846 B CN107292846 B CN 107292846B CN 201710501987 A CN201710501987 A CN 201710501987A CN 107292846 B CN107292846 B CN 107292846B
Authority
CN
China
Prior art keywords
projection
theta
frequency domain
data
equation
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.)
Active
Application number
CN201710501987.9A
Other languages
Chinese (zh)
Other versions
CN107292846A (en
Inventor
徐圆
马健晖
周凌宏
丁振强
黄杰星
罗宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southern Medical University
Original Assignee
Southern Medical University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southern Medical University filed Critical Southern Medical University
Priority to CN201710501987.9A priority Critical patent/CN107292846B/en
Publication of CN107292846A publication Critical patent/CN107292846A/en
Application granted granted Critical
Publication of CN107292846B publication Critical patent/CN107292846B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)

Abstract

The invention relates to a method for recovering incomplete CT projection data under a circular orbit, which firstly deduces an accurate consistency condition under circular orbit scanning data according to a John equation and realizes a corresponding projection data recovery method according to the condition. The method comprises the following steps: firstly, obtaining the position information of a bad pixel point; then, initializing the bad pixel point projection by using a spatial interpolation method; respectively carrying out Fourier transform on the initialized projection and the projections at two adjacent angles in front and at the back of the initialized projection to obtain transformed projection frequency data; carrying out Fourier inversion on the weighted projection frequency domain data to generate corrected projection data, and obtaining a first iteration final value; and taking the final result of the last iteration as an initial value, and repeating the third step to the fifth step until the requirement of the iteration times is met. The invention has the advantages of image quality, simple realization and applicability.

Description

Recovery method of incomplete CT projection data under circular orbit
Technical Field
The invention relates to medical image processing, in particular to image enhancement or restoration, and particularly relates to a method for restoring incomplete CT projection data under a circular orbit.
Background
Cone beam computed tomography is a medical imaging method that can obtain three-dimensional anatomical images of a human body quickly and with low radiation dose. The advent of digital flat panel detectors has made possible the clinical use of true Cone-beam Computed Tomography (CBCT).
Since CBCT employs a cone-beam scanning architecture, the acquired projection data is highly redundant, and Consistency Conditions (Consistency Conditions) must be satisfied to take advantage of the redundancy of the data to achieve high efficiency. John derived a hyper-hyperbolic Partial Differential Equation (PDE) as a consistency condition for line integrals in 1938 from a mathematical point of view, so this Equation is called john's Equation. In the 70 s, with the advent and prosperity of CT, John's Equation was widely used in the CT field and was more deeply derived and extended.
Until now, many scholars have proposed consistency conditions based on the john equation, such as the one proposed by Patch (Sarah k. Patch, calculation of Unmeasured Third-Generation VCT Views From Measured Views,2002, IEEE trans. med. imaging), in which authors have replaced the parameters in the cartesian coordinate system based john equation by using the Third-Generation spiral CT parameters and have performed a series of deep derivations on the basis of the parameters, converted the second order partial differential equation into a set of first order partial differential equations, fourier transformed the equations, and then applied the form of the frequency domain to calculate the Unmeasured projection data From the obtained projections. According to this theory, the literature uses a "rotation + helix + rotation" source trajectory with a pitch of 0.128m, 984 projections are acquired per revolution, and then extrapolated against the acquired projection data using a fourier transform form corresponding to the derivation equation to obtain the unmeasured projection data. However, the consistency condition derived by Patch is derived based on the geometry of spiral CT, and the spiral structure requires that the variable z continuously changes along the longitudinal axis, which introduces a variable into the formula, and makes the formula complex to calculate and difficult to solve. In the literature, due to the scanning mode of "rotation + spiral + rotation", the moving speed of the table needs to be increased and decreased before and after the scanning stage of "spiral", respectively, to ensure the fast transition and smooth transition of data acquisition in the three scanning stages, but in practical application, the moving speed of the table of general spiral CT is constant, and therefore, the consistency condition in the literature is difficult to be applied specifically. In addition, Hao Yan et al used Patch's consistency conditions based on three-generation helical CT when using a moving beam stop array (bead stop array, BSA) for scatter correction (Hao Yan, Xuanqin Mou, Shaojie Tang, Qiong Xu and Maria Zankl, project movingbead stop array,2010, Phys.Med.biol.). The beam stop array is a mask on which the mask spots are regularly arranged. The X-rays are irradiated onto the object through the shielding plate, wherein the X-rays passing through the shielding point on the shielding plate are shielded. If a value is present at a corresponding location of a detector occlusion point, this value is a scatter value. Since the BSA is used, the obtained projection data is incomplete, and thus needs to be restored using a consistency condition. The document adopts a derivation formula [ Patch (2002) ] of Patch based on three generations of spiral CT, and combines a spatial interpolation method to provide a PC _ SI method for the field of scatter correction. However, as mentioned above, since the z-axis of the spiral CT is continuously changed, which makes the formula difficult to be applied, the document refers to the method of Defrise et al, and the z-coordinate is replaced by the vertical coordinate of the detector, but the accuracy is reduced by this approximation. From a geometrical point of view, the spiral CT can be degenerated to a circular track CBCT by only making z equal to 0, but the formula in the literature is to differentiate z, and when z is constant, the derivative has no meaning, so the method used in the literature is not suitable for the circular track CBCT.
Therefore, it is necessary to provide a method for recovering incomplete CT projection data under a circular orbit to overcome the deficiencies of the prior art.
Disclosure of Invention
Compared with the traditional spatial interpolation method, the method can better recover high-frequency components in lost information, further greatly reduce the bar artifacts in the image and obviously improve the image quality. The method of using the Patch under the circular orbit needs to approximate the ordinate of the detector to replace the z coordinate, the consistency condition adopted by the method is obtained by derivation based on the circular orbit CBCT geometric structure, the precision reduction caused by the approximation of the z coordinate can be avoided, and the specific application of the method is difficult because the z axis in the method of the Patch is continuously changed. On the contrary, the method of the invention is only a simple first-order derivation formula, is simple to realize and has applicability.
The invention adopts the following technical scheme:
a method for recovering incomplete CT projection data under a circular orbit is carried out by the following steps:
firstly, obtaining position information of a bad pixel point;
secondly, performing linear interpolation on the bad pixel points in each projection by using a horizontal one-dimensional cubic spline interpolation method to recover low-frequency information of the bad pixel points, and obtaining an initial projection at each angle through the step;
step three, making N equal to 1, entering the step four,
fourthly, enabling M to be equal to N, wherein M is the current iteration number;
performing two-dimensional Fourier transform on the initialized projection at each angle to obtain frequency domain data of the initialized projection at each angle in a frequency domain;
fifthly, respectively bringing the frequency domain forms of the front and the back adjacent projections into a formula corresponding to the consistency condition to calculate a projection result;
sixthly, performing Fourier inverse transformation on the weighted projection to obtain an Mth iteration final value;
seventhly, judging whether the current iteration times M are more than or equal to the iteration termination times S, and if so, entering the ninth step; otherwise, entering the eighth step;
eighthly, taking the M iteration final value as an initialization projection, enabling N to be N +1, and entering the third step;
the ninth step; and taking the final value of the Mth iteration as a finally recovered data result.
Preferably, the above method for recovering incomplete CT projection data under circular orbit,
the consistency condition used in the fifth step is specifically as follows:
Figure BDA0001333901850000031
where ρ is the distance from the source to the center of rotation, d is the distance from the center of rotation to the detector, and (α)1,α2) Is the flat panel detector coordinate and theta is the azimuth angle.
Preferably, in the fifth step, the two previous and next projection results are weighted by using the same weighting factor ω, so as to obtain the projection result after updating the weighting.
Preferably, the above method for recovering incomplete CT projection data under circular orbit,
and fifthly, recovering the frequency domain data of each initialized projection by using a consistency condition, namely a formula (28), wherein the formula (28) is specifically formed as follows:
u*(θ+dθ)＝u*(θ)+dθ·U,k2equation (28) ≠ 0 … …;
Figure BDA0001333901850000032
wherein u (theta + d theta) and u (theta) respectively represent frequency domain data of the initialized projection under theta + d theta and theta angles, u (theta) is a known quantity, u (theta + d theta) is a required unknown quantity, d theta represents the interval between two adjacent angles, k1、k2Is the frequency domain coordinate of the flat panel detector, i represents the imaginary part of the complex number;
the specific method of use of equation (28) is as follows:
u (θ) is first calculated by substituting u (θ + d θ) into equation (28):
u*(θ)＝u*(θ+dθ)-dθ·U；
then, u (θ) in the formula (28) is replaced with u (θ -d θ), and u (θ + d θ) in the formula (28) is replaced with u (θ -d θ) to calculate:
u*(θ)＝u*(θ-dθ)+dθ·U；
as any one of the two adjacent projection frequency domain data before and after selection can realize projection data recovery, in order to balance the recovery effect, a weight factor omega is used, and omega belongs to [0,1 ]]Carrying out weighted average on the calculation results of the two formulas to obtain a weighted result uR*(θ), i.e.:
uR*(θ)＝ω(u*(θ-dθ)+dθ·U)+(1-ω)(u*(θ+dθ)-dθ·U),k2≠0；
since the formula (28) is only at k2Not equal to 0, so that k is valid for the frequency domain data u x (theta)2A portion of 0, directly using the corresponding portion in the initialized projection frequency domain data at the angle θ to compensate;
a sixth step, specifically for uR*(theta) carrying out inverse Fourier transform to obtain the Mth iteration final value uR(θ)；
uR(θ)＝F-1(uR*(θ))……(1)；
Wherein u in formula (1)R(theta) represents frequency domain data u weighted at theta angleR*Result of inverse Fourier transform of (theta), F-1Representing an inverse fourier transform;
step seven, specifically, judging whether the current iteration times M is more than or equal to the iteration termination times S, if so, entering the step ninth; otherwise, entering the eighth step;
the eighth step, specifically, the M-th iteration final value uR(θ) as an initialization projection, making N equal to N +1, and proceeding to a third step;
the ninth step, specifically, the final value u of the Mth iterationR(θ) is the final recovered data result.
Preferably, in the method for recovering incomplete CT projection data under a circular orbit, the first step of obtaining the position information of the bad pixel point is to specifically convert the actually acquired projection data into a corresponding chord graph, that is: the x axis of the chord graph is the position of a projected pixel point, and the y axis is the position of a projection angle corresponding to the pixel position; under the condition that a certain detector unit is damaged, a vertical line appears at the corresponding pixel position in the chord graph, the position information of a bad pixel point is obtained in the chord graph according to the theory, and the position of the bad pixel point is recorded.
The derivation process of the formula (28) applied in the method is as follows:
since the consistency condition proposed by the present invention is based on the john equation, the john equation is briefly introduced and analyzed. FJohn proposed a hyper-hyperbolic second-order partial DIFFERENTIAL EQUATION in 1938 in THE document "THE ultra-HYPERBOLIC DIFFERENTIAL EQUATION WITHURUR INDEPENDENDENT VARIABLES" as a consistency condition of line integrals, and THE specific form is as follows:
Figure BDA0001333901850000041
in the formula, i and j are subscripts;
u is the line integral of the function f across and η:
u(；η)＝∫f(+t(η-))dt (3)；
the following steps are the consistency condition derivation steps based on John's equation under the circular orbit:
firstly, parameterizing a sum eta by adopting a circular orbit CBCT coordinate parameter:
Figure BDA0001333901850000042
Figure BDA0001333901850000043
where ρ is the distance from the source to the center of rotation and d is the distance from the center of rotation to the detector, (α)1,α2) Is the flat panel detector coordinate and theta is the azimuth angle.
Separately summing each component of η to compute the partial derivative:
Figure BDA0001333901850000051
Figure BDA0001333901850000052
Figure BDA0001333901850000053
Figure BDA0001333901850000054
Figure BDA0001333901850000055
Figure BDA0001333901850000056
for writing convenience, the following operators are defined:
Figure BDA0001333901850000057
Figure BDA0001333901850000058
Figure BDA0001333901850000059
Figure BDA00013339018500000510
Figure BDA00013339018500000511
substituting equations (6) to (11) into equation (2) yields:
Figure BDA00013339018500000512
Figure BDA00013339018500000513
Figure BDA00013339018500000514
by the formulae (13) to (16) to H1，H2,
Figure BDA00013339018500000515
Substituting into equations (17) to (19) results in:
Figure BDA0001333901850000061
from the above equation, we derive:
Figure BDA0001333901850000062
merging the same kind of terms for the above formula (21) transposition:
Figure BDA0001333901850000063
Figure BDA0001333901850000064
solving equation (23), the final derivation is as follows:
Figure BDA0001333901850000065
equation (24) is the consistency condition based on John's equation under the circular orbit we derive. The expression form of the corresponding frequency domain obtained by performing fourier transform on the formula (24) is as follows:
Figure BDA0001333901850000066
in the above formula, uθ *Denotes the Fourier transform of u derivative to θ, (k)1,k2) Is (alpha)1,α2) Corresponding frequency domain form of (a).
For ease of expression and application, the right equation of equation (25) is replaced with U:
Figure BDA0001333901850000071
in view of the general definition of the derivation formula, equation (25) can be modified to the form:
Figure BDA0001333901850000072
namely:
u*(θ+dθ)＝u*(θ)+dθ·U,k2≠0 (28)；
equation (28) is an improved consistency condition for the present invention, which uses equation (28) for the recovery of the projection data.
The traditional spatial interpolation technology uses other complete pixel data projected at the same angle to recover lost data, and the used reconstruction information is only limited to the projection of the angle of the lost data, so that the interpolation method can only recover low-frequency information in the lost data, and basically does not recover the lost high-frequency information. The consistency condition applied by the method of the invention is that the lost data is recovered based on the data in the same projection, and the corresponding position data information in the adjacent angle projections is utilized. Therefore, the method can well recover the high-frequency information in the lost data by adopting the consistency condition, can basically eliminate the bar artifact in the reconstructed image obtained by the interpolation technology, and obviously improves the image quality.
Drawings
FIG. 1 is a geometric diagram and a parametric diagram of an exemplary circular CBCT scan, where ρ is the distance from the source to the center of rotation, d is the distance from the center of rotation to the detector, and (α)1,α2) Is the flat panel detector coordinate and theta is the azimuth angle.
Fig. 2 is a chord graph corresponding to projection data, where the x-axis of the abscissa is the pixel position and the y-axis of the ordinate is the projection angle.
FIG. 3 is a flow chart of a data recovery method of the present invention.
Fig. 4 is a schematic diagram of simulation results of three methods, fig. 4(a) is an ideal image, fig. 4(b) is an original image without any processing, fig. 4(c) is a result obtained by linear difference, fig. 4(d), (e), (f) are result diagrams after one iteration, two iterations and three iterations of the method of the present invention, and fig. (g), (h) and (i) are result diagrams after one iteration, two iterations and three iterations of the Patch method.
Detailed Description
Example 1.
A method for recovering incomplete CT projection data under a circular orbit is carried out by the following steps:
the first step, obtain bad pixel point positional information, convert the projection data of actual collection into corresponding chord graph, promptly: the x axis of the chord graph is the position of a projected pixel point, and the y axis is the position of a projection angle corresponding to the pixel position; under the condition that a certain detector unit is damaged, a vertical line appears at the corresponding pixel position in the chord graph, the position information of a bad pixel point is obtained in the chord graph according to the theory, and the position of the bad pixel point is recorded.
And secondly, performing linear interpolation on the bad pixel points in each projection by using a horizontal one-dimensional cubic spline interpolation method to recover the low-frequency information of the bad pixel points, and obtaining the initialized projection at each angle through the step.
Step three, making N equal to 1, entering the step four,
fourthly, enabling M to be equal to N, wherein M is the current iteration number;
performing two-dimensional Fourier transform on the initialized projection at each angle to obtain frequency domain data of the initialized projection at each angle in a frequency domain;
fifthly, respectively bringing the frequency domain forms of the front and the back adjacent projections into a formula corresponding to the consistency condition to calculate a projection result; specifically, the same weighting factor ω is used to weight the two previous and next projection results, so as to obtain the projection result after updating the weighting.
The consistency conditions used in the fifth step are of the form:
Figure BDA0001333901850000081
where ρ is the distance from the source to the center of rotation, d is the distance from the center of rotation to the detector, and (α)1,α2) Is the flat panel detector coordinate and theta is the azimuth angle, as shown in figure 1.
Fifth, the frequency domain data of each initialized projection may be recovered by using a consistency condition, that is, formula (28), where the formula (28) is specifically formed as follows:
u*(θ+dθ)＝u*(θ)+dθ·U,k2equation (28) ≠ 0 … …;
Figure BDA0001333901850000082
wherein u (theta + d theta) and u (theta) respectively represent frequency domain data of the initialized projection under theta + d theta and theta angles, u (theta) is a known quantity, u (theta + d theta) is a required unknown quantity, d theta represents the interval between two adjacent angles, k1、k2Is the frequency domain coordinate of the flat panel detector, i represents the imaginary part of the complex number;
the specific method of use of equation (28) is as follows:
u (θ) is first calculated by substituting u (θ + d θ) into equation (28):
u*(θ)＝u*(θ+dθ)-dθ·U；
then, u (θ) in the formula (28) is replaced with u (θ -d θ), and u (θ + d θ) in the formula (28) is replaced with u (θ -d θ) to calculate:
u*(θ)＝u*(θ-dθ)+dθ·U；
as any one of the two adjacent projection frequency domain data before and after selection can realize projection data recovery, in order to balance the recovery effect, a weight factor omega is used, and omega belongs to [0,1 ]]Carrying out weighted average on the calculation results of the two formulas to obtain a weighted result uR*(θ), i.e.:
uR*(θ)＝ω(u*(θ-dθ)+dθ·U)+(1-ω)(u*(θ+dθ)-dθ·U),k2≠0；
since the formula (28) is only at k2Not equal to 0, so that k is valid for the frequency domain data u x (theta)2A portion of 0, directly using the corresponding portion in the initialized projection frequency domain data at the angle θ to compensate;
a sixth step, specifically for uR*(theta) carrying out inverse Fourier transform to obtain the Mth iteration final value uR(θ)；
uR(θ)＝F-1(uR*(θ))……(1)；
Wherein u in formula (1)R(theta) represents frequency domain data u weighted at theta angleR*Result of inverse Fourier transform of (theta), F-1Representing an inverse fourier transform.
Step seven, specifically, judging whether the current iteration times M is more than or equal to the iteration termination times S, if so, entering the step ninth; otherwise, entering the eighth step.
The eighth step, specifically, the M-th iteration final value uR(θ) as the initialization projection, let N be N +1, and proceed to the third step.
The ninth step, specifically, the final value u of the Mth iterationR(θ) is the final recovered data result.
The derivation process of the formula (28) applied in the method is as follows:
since the consistency condition proposed by the present invention is based on the john equation, the john equation is briefly introduced and analyzed. F John, in 1938, proposed a hyper-hyperbolic second-order partial DIFFERENTIAL EQUATION in THE document "THE ultra-hybrid EQUATION WITH FOUR exponential variations, as a consistency condition for line integrals, in THE following specific form:
Figure BDA0001333901850000091
u is the line integral of the function f across and η:
u(；η)＝∫f(+t(η-))dt (3)；
the following steps are the consistency condition derivation steps based on John's equation under the circular orbit:
firstly, parameterizing a sum eta by adopting a circular orbit CBCT coordinate parameter:
Figure BDA0001333901850000092
Figure BDA0001333901850000093
where ρ is the distance from the source to the center of rotation and d is the distance from the center of rotation to the detector, (α)1,α2) Is the flat panel detector coordinate and theta is the azimuth angle.
Separately summing each component of η to compute the partial derivative:
Figure BDA0001333901850000101
Figure BDA0001333901850000102
Figure BDA0001333901850000103
Figure BDA0001333901850000104
Figure BDA0001333901850000105
Figure BDA0001333901850000106
for writing convenience, the following operators are defined:
Figure BDA0001333901850000107
Figure BDA0001333901850000108
Figure BDA0001333901850000109
Figure BDA00013339018500001010
Figure BDA00013339018500001011
substituting equations (6) to (11) into equation (2) yields:
Figure BDA00013339018500001012
Figure BDA00013339018500001013
Figure BDA00013339018500001014
is represented by the formula (13) to (16)H1，H2,
Figure BDA00013339018500001015
Substituting into equations (17) to (19) results in:
Figure BDA0001333901850000111
from the above equation, we derive:
Figure BDA0001333901850000112
merging the same kind of terms for the above formula (21) transposition:
Figure BDA0001333901850000113
Figure BDA0001333901850000114
solving equation (23), the final derivation is as follows:
Figure BDA0001333901850000115
equation (24) is the consistency condition based on John's equation under the circular orbit we derive. The expression form of the corresponding frequency domain obtained by performing fourier transform on the formula (24) is as follows:
Figure BDA0001333901850000116
in the above formula, uθ *Denotes the Fourier transform of u derivative to θ, (k)1,k2) Is (alpha)1,α2) Corresponding frequency domain form of (a).
For ease of expression and application, the right equation of equation (25) is replaced with U:
Figure BDA0001333901850000121
in view of the general definition of the derivation formula, equation (25) can be modified to the form:
Figure BDA0001333901850000122
namely:
u*(θ+dθ)＝u*(θ)+dθ·U,k2≠0 (28)；
equation (28) is an improved consistency condition for the present invention, which uses equation (28) for the recovery of the projection data.
The traditional spatial interpolation technology uses other complete pixel data projected at the same angle to recover lost data, and the used reconstruction information is only limited to the projection of the angle of the lost data, so that the interpolation method can only recover low-frequency information in the lost data, and basically does not recover the lost high-frequency information. The consistency condition applied by the method of the invention is that the lost data is recovered based on the data in the same projection, and the corresponding position data information in the adjacent angle projections is utilized. Therefore, the method can well recover the high-frequency information in the lost data by adopting the consistency condition, can basically eliminate the bar artifact in the reconstructed image obtained by the interpolation technology, and obviously improves the image quality.
Compared with the traditional spatial interpolation method, the method can well recover the high-frequency components in the lost information, further greatly reduce the bar artifacts in the image and obviously improve the image quality. The method of using the Patch under the circular orbit needs to approximate the ordinate of the detector to replace the z coordinate, the consistency condition adopted by the method is obtained by derivation based on the circular orbit CBCT geometric structure, the precision reduction caused by the approximation of the z coordinate can be avoided, and the specific application of the method is difficult because the z axis in the method of the Patch is continuously changed. On the contrary, the method of the invention is only a simple first-order derivation formula, is simple to realize and has applicability.
Example 2.
In this embodiment, a conventional spatial interpolation method, a Patch method (detailed description in the background art), and the method are respectively used to recover a set of incomplete projection data, where the incomplete projection data is acquired by a CBCT with bad pixel points of a flat panel detector, and CT scanning parameters are as follows: the distance rho from the ray source to the rotation center is 500mm, the distance d from the rotation center to the detector is 500mm, the size of the flat panel detector is 850 multiplied by 200mm, the size of the pixel point is 1 multiplied by 1mm, and 1080 projections are collected around the object in a 360-degree rotation mode.
The specific implementation method of the technology of the invention is as follows:
firstly, obtaining the position information of the bad pixel points. The actually acquired projection data are converted into a corresponding chord graph as shown in fig. 2. The horizontal coordinate corresponding to the black vertical line in the graph is the position of the bad pixel point, and the position information of all the bad pixel points can be found and recorded through the method.
And secondly, performing linear interpolation on the bad pixel points in each projection by using a horizontal one-dimensional cubic spline interpolation method, and taking the interpolated result as the initial value of the bad pixel point.
And thirdly, carrying out Fourier transform on the projection after interpolation and the front and back adjacent projections to obtain a corresponding frequency domain form so as to conveniently apply consistency conditions deduced by the method.
And step four, substituting the frequency domain forms of the two adjacent projections obtained in the previous step into formula (28), recovering the lost information (especially high-frequency component information) of the middle projection, and weighting the two results by using the same weight factor omega, wherein the example makes omega equal to 0.5 because no more constraint conditions indicate which of the two adjacent projections contributes more to the middle projection. The initial projection value is then corrected by the result of the calculation of the weighting.
And fifthly, carrying out inverse Fourier transform on the corrected value to obtain a final value of the first iteration.
And sixthly, taking the final result of the last iteration as an initial value, repeating the third step to the fifth step, and iterating for three times in the embodiment.
The flow chart of the specific steps is shown in fig. 3.
The simulation result is shown in fig. 4, where fig. 4(a) is an ideal image, fig. 4(b) is an original image without any processing, fig. 4(c) is a result obtained by linear difference, fig. 4(d), (e), (f) are result graphs after one iteration, two iterations and three iterations of the method of the present invention, and fig. (g), (h) and (i) are result graphs after one iteration, two iterations and three iterations of the Patch method.
The reconstructed image was analyzed using Absolute Error (Absolute Error):
Figure BDA0001333901850000131
wherein N is the number of pixel points, riIs the gray value corresponding to the ith point of the reconstructed image, tiIs the gray value corresponding to the ith point of the ideal image. The results are as follows:
method of producing a composite material MAE of layer 101
Linear spatial interpolation 1.3368×10+6
One iteration of the Patch method 3.0641×10+5
Patch method iteration 3.1477×10+5
Triple iteration of Patch method 2.6800×10+5
One iteration of the method 2.9724×10+5
Second iteration of the method 2.9658×10+5
Three iterations of the method 2.6187×10+5
From the perspective of the image reconstruction effect, the reconstructed image obtained by using the spatial interpolation method has serious bar artifacts, which seriously affect the extraction of image diagnostic information, and the Patch method and the method can significantly improve the image quality.
Quantitative analysis of the reconstructed image using absolute errors can find: compared with the Patch method, the image quality obtained by the method is closer to the ideal image quality, each iteration effect is better than that of the Patch method, and the later iteration effect is better than that of the previous iteration effect.
Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present invention and not for limiting the protection scope of the present invention, and although the present invention is described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that modifications or equivalent substitutions can be made on the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims (3)

1. A method for recovering incomplete CT projection data under a circular orbit is characterized by comprising the following steps:
firstly, obtaining position information of a bad pixel point;
secondly, performing linear interpolation on the bad pixel points in each projection by using a horizontal one-dimensional cubic spline interpolation method to recover low-frequency information of the bad pixel points, and obtaining an initial projection at each angle through the step;
step three, enabling N to be 1, and entering the step four;
fourthly, enabling M to be equal to N, wherein M is the current iteration number;
performing two-dimensional Fourier transform on the initialized projection at each angle to obtain frequency domain data of the initialized projection at each angle in a frequency domain;
fifthly, respectively bringing the frequency domain forms of the front and the back adjacent projections into a formula corresponding to the consistency condition to calculate a projection result; weighting the front projection result and the rear projection result by using the same weight factor omega to obtain a weighted and updated projection result;
sixthly, performing Fourier inversion on the weighted projection result to obtain an Mth iteration final value;
seventhly, judging whether the current iteration times M are more than or equal to the iteration termination times S, and if so, entering the ninth step; otherwise, entering the eighth step;
eighthly, taking the M iteration final value as an initialization projection, enabling N to be N +1, and entering the third step;
the ninth step, the final value of the Mth iteration is taken as the final recovered data result;
the consistency condition used in the fifth step is specifically as follows:
Figure FDA0002673576870000011
where ρ is the distance from the source to the center of rotation, d is the distance from the center of rotation to the detector, (α)1,α2) Is the flat panel detector coordinate, theta is the azimuth;
u is the line integral of the two points of the target function f crossing and η:
u (; η) ═ f (+ t (η -)) dt; t is the sign of the integral variable at two points where the ray f passes through and η.
2. The method for recovering incomplete CT projection data under a circular orbit as claimed in claim 1, wherein: performing fourier transform on the consistency condition in the fifth step to obtain an expression form of a corresponding frequency domain as follows:
Figure FDA0002673576870000012
wherein u isθ *Denotes the Fourier transform of u derivative to θ, (k)1,k2) Is (alpha)1,α2) The corresponding frequency domain form of (a);
substituting the right equation of equation (25) with U, we get equation (26):
Figure FDA0002673576870000021
wherein u (theta + d theta) and u (theta) respectively represent frequency domain data of the initialized projection under theta + d theta and theta angles, u (theta) is a known quantity, u (theta + d theta) is a required unknown quantity, d theta represents the interval between two adjacent angles, and (k) represents the interval between two adjacent angles1,k2) Is the frequency domain coordinate of the flat panel detector, i represents the imaginary part of the complex number;
rewrite equation (25) to the following form:
Figure FDA0002673576870000022
equation deformation of equation (27) yields the following form:
u*(θ+dθ)＝u*(θ)+dθ·U,k2≠0 (28)；
fifthly, restoring the frequency domain data of each initialized projection by using a formula (28);
the specific method of use of equation (28) is as follows:
u (θ) is first calculated by substituting u (θ + d θ) into equation (28):
u*(θ)＝u*(θ+dθ)-dθ·U；
then, u (θ) in the formula (28) is replaced with u (θ -d θ), and u (θ + d θ) in the formula (28) is replaced with u (θ -d θ) to calculate:
u*(θ)＝u*(θ-dθ)+dθ·U；
as any one of the two adjacent projection frequency domain data before and after selection can realize projection data recovery, in order to balance the recovery effect, a weight factor omega is used, and omega belongs to [0,1 ]]Carrying out weighted average on the calculation results of the two formulas to obtain a weighted result uR*(θ), i.e.:
uR*(θ)＝ω(u*(θ-dθ)+dθ·U)+(1-ω)(u*(θ+dθ)-dθ·U),k2≠0；
since the formula (28) is only at k2Not equal to 0, so that k is valid for the frequency domain data u x (theta)2A portion of 0, directly using the corresponding portion in the initialized projection frequency domain data at the angle θ to compensate;
a sixth step, specifically for uR*(theta) carrying out inverse Fourier transform to obtain the Mth iteration final value uR(θ)；
uR(θ)＝F-1(uR*(θ)) (1)；
Wherein u in formula (1)R(theta) represents frequency domain data u weighted at theta angleR*Result of inverse Fourier transform of (theta), F-1Representing an inverse fourier transform;
step seven, specifically, judging whether the current iteration times M is more than or equal to the iteration termination times S, if so, entering the step ninth; otherwise, entering the eighth step;
the eighth step, specifically, the M-th iteration final value uR(θ) as an initialization projection, making N equal to N +1, and proceeding to a third step;
the ninth step, specifically, the final value u of the Mth iterationR(θ) is the final recovered data result.
3. The method of claim 2, wherein the incomplete CT projection data under the circular orbit is recovered,
the first step, obtaining the position information of the bad pixel point is to convert the actually collected projection data into a corresponding chord graph, namely: the x axis of the chord graph is the position of a projected pixel point, and the y axis is the position of a projection angle corresponding to the pixel position; under the condition that a certain detector unit is damaged, a vertical line appears at the corresponding pixel position in the chord graph, the position information of a bad pixel point is obtained in the chord graph according to the theory, and the position of the bad pixel point is recorded.
CN201710501987.9A 2017-06-27 2017-06-27 Recovery method of incomplete CT projection data under circular orbit Active CN107292846B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710501987.9A CN107292846B (en) 2017-06-27 2017-06-27 Recovery method of incomplete CT projection data under circular orbit

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710501987.9A CN107292846B (en) 2017-06-27 2017-06-27 Recovery method of incomplete CT projection data under circular orbit

Publications (2)

Publication Number Publication Date
CN107292846A CN107292846A (en) 2017-10-24
CN107292846B true CN107292846B (en) 2020-11-10

Family

ID=60098750

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710501987.9A Active CN107292846B (en) 2017-06-27 2017-06-27 Recovery method of incomplete CT projection data under circular orbit

Country Status (1)

Country Link
CN (1) CN107292846B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6173030B1 (en) * 1998-11-25 2001-01-09 General Electric Company Almost-everywhere extrapolation using 2D transforms from cone-beam data
EP1769460B1 (en) * 2004-07-13 2008-01-16 Philips Intellectual Property & Standards GmbH A computed tomography method for the reconstruction of object images from real and fictitious measured values
CN101627919A (en) * 2009-08-20 2010-01-20 浙江大学 PET concentration reestablishing method based on Kalman filtration in limited sampling angle
CN102680948A (en) * 2012-05-15 2012-09-19 东南大学 Method for estimating modulation frequency and starting frequency of linear frequency-modulated signal
CN104021582A (en) * 2014-05-28 2014-09-03 山东大学 CT (Computed Tomography) iterative image reconstruction method
CN104050631A (en) * 2013-11-25 2014-09-17 中国科学院上海应用物理研究所 Low-dose CT image reconstruction method
CN105590332A (en) * 2015-12-24 2016-05-18 电子科技大学 Rapid algebraic reconstruction technique applied to computed tomography imaging
CN106612427A (en) * 2016-12-29 2017-05-03 浙江工商大学 Method for generating spatial-temporal consistency depth map sequence based on convolution neural network

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6173030B1 (en) * 1998-11-25 2001-01-09 General Electric Company Almost-everywhere extrapolation using 2D transforms from cone-beam data
EP1769460B1 (en) * 2004-07-13 2008-01-16 Philips Intellectual Property & Standards GmbH A computed tomography method for the reconstruction of object images from real and fictitious measured values
CN101627919A (en) * 2009-08-20 2010-01-20 浙江大学 PET concentration reestablishing method based on Kalman filtration in limited sampling angle
CN102680948A (en) * 2012-05-15 2012-09-19 东南大学 Method for estimating modulation frequency and starting frequency of linear frequency-modulated signal
CN104050631A (en) * 2013-11-25 2014-09-17 中国科学院上海应用物理研究所 Low-dose CT image reconstruction method
CN104021582A (en) * 2014-05-28 2014-09-03 山东大学 CT (Computed Tomography) iterative image reconstruction method
CN105590332A (en) * 2015-12-24 2016-05-18 电子科技大学 Rapid algebraic reconstruction technique applied to computed tomography imaging
CN106612427A (en) * 2016-12-29 2017-05-03 浙江工商大学 Method for generating spatial-temporal consistency depth map sequence based on convolution neural network

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Consistency conditions upon 3D CT data and the wave equation;S K Patch 等;《Physics in Medicine and Biology》;20021231(第47期);第1-14页 *
Improved 2D rebinning of helical cone-beam CT data using John"s equation;M. Defrise 等;《2002 IEEE Nuclear Science Symposium Conference Record》;20031027;第1465-1469页 *
WE-AB-207A-02:John’s Equation Based Consistency Condition and Incomplete Projection Restoration Upon Circular Orbit CBCT;J. Ma 等;《Medical Physics》;20160630;第43卷(第6期);参见第3797-3798页WE-AB-207A-02部分 *
基于投影域数据恢复的低剂量CT稀疏角度重建;刘文磊 等;《CT理论与应用研究》;20130731;第22卷(第3期);第421-428页 *

Also Published As

Publication number Publication date
CN107292846A (en) 2017-10-24

Similar Documents

Publication Publication Date Title
Bao et al. Convolutional sparse coding for compressed sensing CT reconstruction
Xu et al. A practical cone-beam CT scatter correction method with optimized Monte Carlo simulations for image-guided radiation therapy
Gjesteby et al. A dual-stream deep convolutional network for reducing metal streak artifacts in CT images
Niu et al. Sparse-view x-ray CT reconstruction via total generalized variation regularization
Wang et al. FBP-Net for direct reconstruction of dynamic PET images
CN109146994B (en) Metal artifact correction method for multi-energy spectrum X-ray CT imaging
Yan et al. Projection correlation based view interpolation for cone beam CT: primary fluence restoration in scatter measurement with a moving beam stop array
CN111292386B (en) CT projection metal trace complement metal artifact correction method based on U-net
CN104700438B (en) Image rebuilding method and device
CN103810734B (en) A kind of low dose X-ray CT data for projection restoration methods
CN104408758A (en) Low-dose processing method of energy spectrum CT image
CN103479379B (en) A kind of image rebuilding method of tilting screw scanning and device
Wu et al. Iterative CT shading correction with no prior information
Chen et al. 4D-AirNet: a temporally-resolved CBCT slice reconstruction method synergizing analytical and iterative method with deep learning
CN104751429A (en) Dictionary learning based low-dosage energy spectrum CT image processing method
Zhang et al. Sparse-view X-ray CT reconstruction with Gamma regularization
Yang et al. Shading correction assisted iterative cone-beam CT reconstruction
CN116452423A (en) Simultaneous sparse angle CT reconstruction and metal artifact high-precision correction method
CN103793890A (en) Method for recovering and processing energy spectrum CT images
Zhang et al. 4D radiomics: impact of 4D-CBCT image quality on radiomic analysis
Biguri et al. A general method for motion compensation in x-ray computed tomography
Zhang et al. Dynamic cone-beam CT reconstruction using spatial and temporal implicit neural representation learning (STINR)
CN105844678A (en) Low dose X-ray CT image reconstruction method based on completely generalized variational regularization
Qiu et al. New iterative cone beam CT reconstruction software: parameter optimisation and convergence study
Kole et al. Evaluation of the ordered subset convex algorithm for cone-beam CT

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