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

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

A 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 in particular to a method for restoring incomplete CT projection data under circular orbits.

Background technique

Cone beam computed tomography is a medical imaging method that can obtain three-dimensional anatomical images of the human body quickly and with low radiation dose. The advent of digital flat detectors makes it possible for true cone beam computed tomography (Cone-beamComputed Tomography, CBCT) to be applied in clinic.

Since CBCT adopts a cone beam scanning structure, the acquired projection data is highly redundant. Consistency Conditions must be satisfied if the data redundancy is to be used to achieve high efficiency. In 1938, John deduced a hyperbolic partial differential equation (Ultrahyperbolic PartialDifferential Equation, PDE) from a mathematical point of view as the consistency condition of the line integral, so the equation is called John's equation. In the 1970s, with the emergence and rise of CT, John's Equation was widely used in the field of CT and was further derived and extended.

Up to now, many scholars have proposed a variety of consistency conditions based on John's equation, such as Patch proposed a consistency condition based on the third generation spiral CT (Sarah K. Views, 2002, IEEE Trans.Med.Imaging), in this paper, the author replaces the parameters in the John's equation based on the Cartesian coordinate system by using the third-generation helical CT parameters, and based on this, a series of In-depth derivation, converting a second-order partial differential equation into a set of first-order ordinary differential equations, Fourier transforming this equation, and then applying its frequency-domain form can be obtained from the projections obtained Unmeasured projection data is calculated. According to this theory, the "rotation + helix + rotation" ray source orbit is used in the literature, the pitch is 0.128m, 984 projections are collected per revolution, and then the obtained projection data is extrapolated using the Fourier transform corresponding to the derivation equation. , to obtain unmeasured projection data. However, the consistency condition derived by Patch is derived based on the geometric structure of the helical CT. The helical structure requires the variable z to change continuously along the vertical axis, which will introduce a change in the formula, which makes the calculation of the formula complicated and solves the problem. Problem is difficult. In the literature, due to the "rotation + helix + rotation" scanning method, it is necessary to increase and decrease the moving speed of the table before and after the "spiral" scanning stage to ensure the fast conversion and smoothness of data acquisition in these three scanning stages. However, in practical applications, the moving speed of the working table of general helical CT is constant, so the consistency conditions in the literature are difficult to be applied concretely. In addition, Hao Yan et al. adopted the consistency condition of Patch based on three generations of helical CT when using a moving beam stop array (BSA) for scatter correction (Hao Yan, Xuanqin Mou, Shaojie Tang, Qiong Xu and Maria Zankl). , Projection movingbeam stop array, 2010, Phys.Med.Biol.). A beam stop array is a shielding plate with regularly arranged shielding points. The X-rays are irradiated onto the object through the shielding plate, and the X-rays passing through the shielding points on the shielding plate will be shielded. If a value appears at the corresponding position of the detector occlusion point, this value is the scattering value. Due to the use of BSA, the resulting projection data is incomplete and needs to be recovered using a consistency condition. This document adopts Patch's derivation formula based on three-generation helical CT [Patch (2002)], combined with the spatial interpolation method, and proposes a PC_SI method for the field of scatter correction. However, as mentioned above, since the z-axis of helical CT is constantly changing, this formula is difficult to apply specifically, so this document refers to the practice of Defrise et al., and uses the ordinate of the detector to approximately replace the z-coordinate. However, this approximation will lead to a decrease in accuracy. From the geometrical point of view, it is only necessary to set z=0, and the spiral CT can be degenerated into a circular orbit CBCT, but the formula in the literature is to derive the z derivative. When z is a constant, the derivative is meaningless, so the method used in this document Not suitable for circular orbit CBCT.

Therefore, in view of the deficiencies of the prior art, it is necessary to provide a method for restoring incomplete CT projection data under a circular orbit to overcome the deficiencies of the prior art.

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is to provide a method for restoring incomplete CT projection data under a circular orbit. Compared with the traditional spatial interpolation method, the method of the present invention can restore the high-frequency components in the lost information well, and further Dramatically reduces striping artifacts in images, significantly improving image quality. The method of using Patch in the circular orbit needs to use the ordinate of the detector to approximate the z coordinate first, and the consistency condition adopted by this method is derived based on the CBCT geometry of the circular orbit, which can avoid the approximation caused by the z coordinate. The accuracy decreases, and the z-axis in the Patch method is constantly changing, making the specific application of this method difficult. On the contrary, the method of the present invention is only a simple first-order derivation formula, which is simple to implement and has applicability.

The present invention adopts following technical scheme:

A method for recovering incomplete CT projection data under a circular orbit is carried out through the following steps:

The first step is to obtain the location information of bad pixels;

In the second step, the horizontal one-dimensional cubic spline interpolation method is used to linearly interpolate the bad pixels in each projection to restore the low-frequency information of the bad pixels. After this step, the initial projection at each angle is obtained;

The third step, let N=1, enter the fourth step,

The fourth step, let M=N, M is the current iteration number;

Perform a two-dimensional Fourier transform on the initialized projection at each angle to obtain the frequency domain data of the initialized projection at each angle in the frequency domain;

In the fifth step, the frequency domain forms of the two adjacent projections before and after are respectively brought into the formula corresponding to the consistency condition to calculate the projection result;

The sixth step is to perform inverse Fourier transform on the weighted projection to obtain the final value of the Mth iteration;

The seventh step is to judge whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step;

The eighth step, take the final value of the Mth iteration as the initialization projection, let N=N+1, and enter the third step;

The ninth step: take the final value of the M-th iteration as the final recovered data result.

Preferably, the recovery method of incomplete CT projection data under the above-mentioned circular orbit,

The consistency condition used in the fifth step is as follows:

where ρ is the distance from the ray source to the rotation center, d is the distance from the rotation center to the detector, (α ₁ , α ₂ ) is the flat panel detector coordinate, and θ is the azimuth angle.

Preferably, the fifth step above specifically uses the same weighting factor ω to weight the two projection results before and after, to obtain a weighted and updated projection result.

Preferably, the recovery method of incomplete CT projection data under the above-mentioned circular orbit,

The fifth step is to restore the frequency domain data of each initialized projection using the consistency condition, that is, formula (28). The specific form of formula (28) is as follows:

u ^* (θ+dθ)=u ^* (θ)+dθ·U,k ₂ ≠0…Formula (28);

Among them, u*(θ+dθ) and u*(θ) represent the frequency domain data of the initialized projection at θ+dθ and θ angles, respectively, u*(θ) is a known quantity, and u*(θ+dθ) is the requirement The unknown quantity of , dθ represents the interval between two adjacent angles, k ₁ and k ₂ are the frequency domain coordinates of the flat panel detector, and i represents the imaginary part of the complex number;

The specific usage of formula (28) is as follows:

First, substitute u*(θ+dθ) into formula (28) to calculate u*(θ):

u ^* (θ)=u ^* (θ+dθ)-dθ·U;

Then use u*(θ-dθ) to replace u*(θ) in formula (28), and u*(θ) to replace u*(θ+dθ) in formula (28) to calculate:

u ^* (θ)=u ^* (θ-dθ)+dθ·U;

Since the projection data recovery can be achieved by selecting any one of the two adjacent projection frequency domain data before and after, in order to balance the recovery effect, the weighting factor ω, ω∈[0,1] is used to perform a weighted average of the calculation results of the above two formulas , the weighted result u ^R* (θ) is obtained, namely:

u ^R* (θ)=ω(u ^* (θ-dθ)+dθ·U)+(1-ω)(u ^* (θ+dθ)-dθ·U),k ₂ ≠0;

Since formula (28) is only valid when k ₂ ≠0, for the part where k ₂ =0 in the frequency domain data u*(θ), directly use the corresponding part in the initialized projected frequency domain data under the angle θ to compensate;

The sixth step is to perform inverse Fourier transform on u ^R* (θ) to obtain the final value u ^R (θ) of the Mth iteration;

u ^R (θ)=F ^-1 (u ^R* (θ))...(1);

Wherein u ^R (θ) in formula (1) represents the inverse Fourier transform result of the weighted frequency domain data u ^R* (θ) under the angle θ, and F ^-1 represents the inverse Fourier transform;

The seventh step is to determine whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step;

The eighth step is to use the final value u ^R (θ) of the Mth iteration as the initialization projection, let N=N+1, and enter the third step;

The ninth step, specifically, takes the final value u ^R (θ) of the Mth iteration as the final recovered data result.

Preferably, in the above-mentioned method for recovering incomplete CT projection data under a circular orbit, the first step is to obtain the position information of bad pixels, specifically converting the actually collected projection data into a corresponding chord diagram, that is, the x-axis of the chord diagram is The pixel position of the projection, the y-axis is the projection angle position corresponding to the pixel position; under the condition that a certain detector unit is damaged, a vertical line will appear on the corresponding pixel position in the chord diagram. According to this theory, in The position information of bad pixels is obtained from the chord diagram, and the positions of bad pixels are recorded.

The derivation process of formula (28) applied in this 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 first. In 1938, FJohn proposed a hyperbolic second-order partial differential equation in the document "THE ULTRAHYPERBOLIC DIFFERENTIAL EQUATION WITHFOURINDEPENDENT VARIABLES" as the consistency condition of the line integral. The specific form is as follows:

In the formula, i and j are subscripts;

u is the line integral of the function f through ε and η:

u(ε; η)=∫f(ε+t(η-ε))dt (3);

The following are the derivation steps of the consistency condition based on the John equation under the circular orbit:

First, ε and η are parameterized using the circular orbit CBCT coordinate parameters:

where ρ is the distance from the ray source to the center of rotation, d is the distance from the center of rotation to the detector, (α ₁ , α ₂ ) are the coordinates of the flat panel detector, and θ is the azimuth angle.

Find the partial derivatives for each component of ε and η separately:

For the convenience of writing, the following operators are defined:

Substitute formulas (6)-(11) into formula (2) to get:

的定义，代入到公式(17)-(19)中，得到：By formulas (13)-(16) for H ₁ , H ₂ ,

The definition of , and substituting it into formulas (17)-(19), we get:

From the above formula, we get:

For the above formula (21), shift items and merge similar items:

Solving formula (23), the final derivation result is as follows:

Equation (24) is the consistency condition based on the John equation under the circular orbit we deduce. Fourier transform is performed on formula (24), and the corresponding expression in the frequency domain is obtained as follows:

In the above formula, u _θ ^* represents the Fourier transform of u after derivation of θ, and (k ₁ , k ₂ ) is the corresponding frequency domain form of (α ₁ , α ₂ ).

For the convenience of expression and application, the right-hand side of formula (25) is replaced by U:

Considering the general definition of the derivation formula, formula (25) can be changed into the following form:

which is:

u ^* (θ+dθ)=u ^* (θ)+dθ·U,k ₂ ≠0 (28);

The formula (28) is the improved consistency condition of the present invention, and the present invention uses the formula (28) to recover the projection data.

The traditional spatial interpolation technology uses other complete pixel data projected at the same angle to restore the lost data, and the reconstruction information used is limited to the projection at the angle of the lost data, so the interpolation method can only restore the low-frequency information in the lost data, There is basically no recovery of the lost high-frequency information. The consistency condition applied by the method of the present invention is not only to recover the lost data based on the data in the same projection, but also to use the corresponding position data information in the adjacent angle projections. Therefore, by adopting the consistency condition, the present invention can restore the high-frequency information in the lost data well, basically eliminate the stripe artifacts in the reconstructed image obtained by using the interpolation technology, and significantly improve the image quality.

Description of drawings

Figure 1 is a schematic diagram of the geometry and parameters of the circular orbit CBCT scanning in the example. In the figure, ρ is the distance from the ray source to the rotation center, d is the distance from the rotation center to the detector, and (α ₁ , α ₂ ) are the coordinates of the flat panel detector. , θ is the azimuth angle.

Figure 2 is the chord diagram corresponding to the projection data, 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 the data recovery method of the present invention.

Figure 4 is a schematic diagram of the simulation results of the three methods, Figure 4(a) is the ideal image, Figure 4(b) is the original image without any processing, Figure 4(c) is the result obtained by linear difference, Figure 4( d), (e), (f) are the results of the first iteration, the second iteration, and the third iteration of the method of the present invention. Figures (g), (h), and (i) are the first iteration, the second iteration, and the third iteration of the Patch method. Result plot after three iterations.

Detailed ways

Example 1.

A method for recovering incomplete CT projection data under a circular orbit is carried out through the following steps:

The first step is to obtain the bad pixel position information, and convert the actually collected projection data into the corresponding chord diagram, that is: the x-axis of the chord diagram is the projected pixel position, and the y-axis is the projection angle position corresponding to the pixel position; Under the condition of damage to a certain detector unit, a vertical line will appear on the corresponding pixel position in the chord diagram. According to this theory, the position information of the bad pixel point is obtained in the chord diagram, and the position of the bad pixel point is recorded.

In the second step, the horizontal one-dimensional cubic spline interpolation method is used to linearly interpolate the bad pixels in each projection to restore the low-frequency information of the bad pixels. After this step, the initial projection at each angle is obtained.

The third step, let N=1, enter the fourth step,

The fourth step, let M=N, M is the current iteration number;

Perform a two-dimensional Fourier transform on the initialized projection at each angle to obtain the frequency domain data of the initialized projection at each angle in the frequency domain;

In the fifth step, the frequency domain forms of the two adjacent projections before and after are respectively brought into the formula corresponding to the consistency condition to calculate the projection result; specifically, the same weight factor ω is used to weight the two projection results before and after, and the weighted update is obtained. The projection result after.

The consistency condition used in the fifth step is of the form:

where ρ is the distance from the ray source to the center of rotation, d is the distance from the center of rotation to the detector, (α ₁ , α ₂ ) is the coordinates of the flat panel detector, and θ is the azimuth angle, as shown in Figure 1.

In the fifth step, the consistency condition, that is, formula (28), can be used to restore the frequency domain data of each initialized projection. The specific form of formula (28) is as follows:

u ^* (θ+dθ)=u ^* (θ)+dθ·U,k ₂ ≠0…Formula (28);

Among them, u*(θ+dθ) and u*(θ) represent the frequency domain data of the initialized projection at θ+dθ and θ angles, respectively, u*(θ) is a known quantity, and u*(θ+dθ) is the requirement The unknown quantity of , dθ represents the interval between two adjacent angles, k ₁ and k ₂ are the frequency domain coordinates of the flat panel detector, and i represents the imaginary part of the complex number;

The specific usage of formula (28) is as follows:

First, substitute u*(θ+dθ) into formula (28) to calculate u*(θ):

u ^* (θ)=u ^* (θ+dθ)-dθ·U;

Then use u*(θ-dθ) to replace u*(θ) in formula (28), and u*(θ) to replace u*(θ+dθ) in formula (28) to calculate:

u ^* (θ)=u ^* (θ-dθ)+dθ·U;

Since the projection data recovery can be achieved by selecting any one of the two adjacent projection frequency domain data before and after, in order to balance the recovery effect, the weighting factor ω, ω∈[0,1] is used to perform a weighted average of the calculation results of the above two formulas , the weighted result u ^R* (θ) is obtained, namely:

u ^R* (θ)=ω(u ^* (θ-dθ)+dθ·U)+(1-ω)(u ^* (θ+dθ)-dθ·U),k ₂ ≠0;

Since formula (28) is only valid when k ₂ ≠0, for the part where k ₂ =0 in the frequency domain data u*(θ), directly use the corresponding part in the initialized projected frequency domain data under the angle θ to compensate;

The sixth step is to perform inverse Fourier transform on u ^R* (θ) to obtain the final value u ^R (θ) of the Mth iteration;

u ^R (θ)=F ^-1 (u ^R* (θ))...(1);

In the formula (1), u ^R (θ) represents the inverse Fourier transform result of the weighted frequency domain data u ^R* (θ) at the angle θ, and F ⁻¹ represents the inverse Fourier transform.

The seventh step is to determine whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step.

The eighth step is to take the final value u ^R (θ) of the Mth iteration as the initialization projection, and let N=N+1, and enter the third step.

The ninth step, specifically, takes the final value u ^R (θ) of the Mth iteration as the final recovered data result.

The derivation process of formula (28) applied in this 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 first. In 1938, F John proposed a hyperbolic second-order partial differential equation in the document "THE ULTRAHYPERBOLIC DIFFERENTIAL EQUATION WITHFOUR INDEPENDENT VARIABLES" as the consistency condition of the line integral. The specific form is as follows:

u is the line integral of the function f through ε and η:

u(ε; η)=∫f(ε+t(η-ε))dt (3);

The following are the derivation steps of the consistency condition based on the John equation under the circular orbit:

First, ε and η are parameterized using the circular orbit CBCT coordinate parameters:

where ρ is the distance from the ray source to the center of rotation, d is the distance from the center of rotation to the detector, (α ₁ , α ₂ ) are the coordinates of the flat panel detector, and θ is the azimuth angle.

Find the partial derivatives for each component of ε and η separately:

For the convenience of writing, the following operators are defined:

Substitute formulas (6)-(11) into formula (2) to get:

的定义，代入到公式(17)-(19)中，得到：By formulas (13)-(16) for H ₁ , H ₂ ,

The definition of , and substituting it into formulas (17)-(19), we get:

From the above formula, we get:

For the above formula (21), shift items and merge similar items:

Solving formula (23), the final derivation result is as follows:

Equation (24) is the consistency condition based on the John equation under the circular orbit we deduce. Fourier transform is performed on formula (24), and the corresponding expression in the frequency domain is obtained as follows:

In the above formula, u _θ ^* represents the Fourier transform of u after derivation of θ, and (k ₁ , k ₂ ) is the corresponding frequency domain form of (α ₁ , α ₂ ).

For the convenience of expression and application, the right-hand side of formula (25) is replaced by U:

Considering the general definition of the derivation formula, formula (25) can be changed into the following form:

which is:

u ^* (θ+dθ)=u ^* (θ)+dθ·U,k ₂ ≠0 (28);

The formula (28) is the improved consistency condition of the present invention, and the present invention uses the formula (28) to recover the projection data.

The traditional spatial interpolation technology uses other complete pixel data projected at the same angle to restore the lost data, and the reconstruction information used is limited to the projection at the angle of the lost data, so the interpolation method can only restore the low-frequency information in the lost data, There is basically no recovery of the lost high-frequency information. The consistency condition applied by the method of the present invention is not only to recover the lost data based on the data in the same projection, but also to use the corresponding position data information in the adjacent angle projections. Therefore, by adopting the consistency condition, the present invention can restore the high-frequency information in the lost data well, basically eliminate the stripe artifacts in the reconstructed image obtained by using the interpolation technology, and significantly improve the image quality.

Compared with the traditional spatial interpolation method, the method of the present invention can restore the high-frequency components in the lost information well, thereby greatly reducing the stripe artifacts in the image and significantly improving the image quality. The method of using Patch in the circular orbit needs to use the ordinate of the detector to approximate the z coordinate first, and the consistency condition adopted by this method is derived based on the CBCT geometry of the circular orbit, which can avoid the approximation caused by the z coordinate. The accuracy decreases, and the z-axis in the Patch method is constantly changing, making the specific application of this method difficult. On the contrary, the method of the present invention is only a simple first-order derivation formula, which is simple to implement and has applicability.

Example 2.

This example uses the traditional spatial interpolation method, the Patch method (detailed in the background art), and the present method to recover a set of incomplete projection data. The incomplete projection data is generated by a flat panel detector with bad pixels. The CT scanning parameters are as follows: the distance ρ 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×200mm, and the pixel size is 1×1mm , rotate 360 degrees around the object to collect 1080 projections.

The specific embodiment of the method of the present technology is as follows:

The first step is to obtain the location information of bad pixels. Convert the actual collected projection data into the corresponding chord diagram, as shown in Figure 2. The abscissa corresponding to the black vertical line in the figure is the position of the bad pixel. In this way, the position information of all the bad pixels can be found and recorded.

In the second step, use the horizontal one-dimensional cubic spline interpolation method to linearly interpolate the bad pixels in each projection, and use the interpolation result as the initial value of the bad pixel.

In the third step, Fourier transform is performed on the interpolated projection and the two adjacent projections before and after to obtain the corresponding frequency domain form, so as to facilitate the application of the consistency condition derived by the present invention.

The fourth step is to substitute the frequency domain form of the two adjacent projections obtained in the previous step into formula (28) to restore the lost information (especially the high-frequency component information) of the intermediate projection, and use the same weight factor ω to The two results are weighted, and since there are no more constraints to indicate which of the two adjacent projections before and after has a greater contribution to the middle projection, this example sets ω = 0.5. The projected initial value is then corrected with the calculated weighted result.

The fifth step is to perform inverse Fourier transform on the corrected value to obtain the final value of the first iteration.

In the sixth step, the final result of the previous iteration is used as the initial value, and the third to fifth steps above are repeated. This example is iterated three times in total.

The flow chart of the specific steps is shown in Figure 3.

The simulation results are shown in Fig. 4. Fig. 4(a) is the ideal image, Fig. 4(b) is the original image without any processing, Fig. 4(c) is the result obtained by linear difference, Fig. 4(d), (e) and (f) are the results of the first iteration, the second iteration and the third iteration of the method of the present invention. Figures (g), (h) and (i) are the first iteration, the second iteration and the third iteration of the Patch method. result graph.

Analyze the reconstructed image using Absolute Error:

where N is the number of pixels, ri is the gray value corresponding to the _ith point of the reconstructed image, and t _i is the gray value corresponding to the ith point of the ideal image. The result is as follows:

method MAE at Tier 101 Linear Space Interpolation 1.3368×10⁺⁶ Patch method one iteration 3.0641×10⁺⁵ Patch method second iteration 3.1477×10⁺⁵ Patch method three iterations 2.6800×10⁺⁵ One iteration of this method 2.9724×10⁺⁵ Second iteration of this method 2.9658×10⁺⁵ This method iterates three times 2.6187×10⁺⁵

From the perspective of the reconstructed image effect, the reconstructed image obtained by using the spatial interpolation method has serious stripe artifacts, which seriously affects the extraction of image diagnostic information, while the Patch method and this method can significantly improve the image quality.

Using the absolute error to quantitatively analyze the reconstructed image, it can be found that: compared with the Patch method, the image quality obtained by the method of the present invention is closer to the ideal image quality, and the effect of each iteration is better than that of the Patch method, and the latter iteration effect is better than the previous one. One iteration effect.

Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit the protection scope of the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should The technical solutions of the present invention may be modified or equivalently replaced without departing from the spirit and scope of the technical solutions of the present invention.

Claims (3)

1. the recovery method of incomplete CT projection data under a circular track, is characterized in that, is carried out by the following steps: The first step is to obtain the location information of bad pixels; In the second step, the horizontal one-dimensional cubic spline interpolation method is used to linearly interpolate the bad pixels in each projection to restore the low-frequency information of the bad pixels. After this step, the initial projection at each angle is obtained; The third step, let N=1, enter the fourth step; The fourth step, let M=N, M is the current number of iterations; Perform a two-dimensional Fourier transform on the initialized projection at each angle to obtain the frequency domain data of the initialized projection at each angle in the frequency domain; In the fifth step, the frequency domain forms of the two adjacent projections before and after are respectively brought into the formula corresponding to the consistency condition to calculate the projection result; the same weighting factor ω is used to weight the two projection results before and after, and the weighted update is obtained. projection result; The sixth step is to perform inverse Fourier transform on the weighted projection result to obtain the final value of the Mth iteration; The seventh step is to judge whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step; The eighth step, take the final value of the Mth iteration as the initialization projection, let N=N+1, and enter the third step; The ninth step, with the final value of the M-th iteration as the final recovered data result; The consistency condition used in the fifth step is as follows:

Among them, ρ is the distance from the ray source to the rotation center, d is the distance from the rotation center to the detector, (α ₁ , α ₂ ) are the coordinates of the flat panel detector, and θ is the azimuth angle; u is the line integral of the objective function f through two points ε and η: u(ε; η)=∫f(ε+t(η-ε))dt; t is the integral variable sign of the ray f passing through the two points ε and η. 2. the recovery method of incomplete CT projection data under the circular orbit according to claim 1, is characterized in that: carry out Fourier transform to the consistency condition in the 5th step, obtain the expression form of corresponding frequency domain as follows:

Among them, u _θ ^* represents the Fourier transform of u after derivation of θ, (k ₁ , k ₂ ) is the corresponding frequency domain form of (α ₁ , α ₂ ); Substitute the right-hand side of formula (25) with U to obtain formula (26):

Among them, u*(θ+dθ) and u*(θ) represent the frequency domain data of the initialized projection at θ+dθ and θ angles, respectively, u*(θ) is a known quantity, and u*(θ+dθ) is the requirement The unknown quantity of , dθ represents the interval between two adjacent angles, (k ₁ , k ₂ ) is the frequency domain coordinate of the flat panel detector, i represents the imaginary part of the complex number; Rewrite formula (25) into the following form:

Equationally transforming formula (27), the following form is obtained: u ^* (θ+dθ)=u ^* (θ)+dθ·U,k ₂ ≠0 (28); The 5th step, specifically uses formula (28) to restore the frequency domain data of each initialization projection; The specific usage of formula (28) is as follows: First, substitute u*(θ+dθ) into formula (28) to calculate u*(θ): u ^* (θ)=u ^* (θ+dθ)-dθ·U; Then use u*(θ-dθ) to replace u*(θ) in formula (28), and u*(θ) to replace u*(θ+dθ) in formula (28) to calculate: u ^* (θ)=u ^* (θ-dθ)+dθ·U; Since the projection data recovery can be achieved by selecting any one of the two adjacent projection frequency domain data before and after, in order to balance the recovery effect, the weighting factor ω, ω∈[0,1] is used to perform a weighted average of the calculation results of the above two formulas , the weighted result u ^R* (θ) is obtained, namely: u ^R* (θ)=ω(u ^* (θ-dθ)+dθ·U)+(1-ω)(u ^* (θ+dθ)-dθ·U),k ₂ ≠0; Since formula (28) is only valid when k ₂ ≠0, for the part where k ₂ =0 in the frequency domain data u*(θ), directly use the corresponding part in the initialized projected frequency domain data under the angle θ to compensate; The sixth step is to perform inverse Fourier transform on u ^R* (θ) to obtain the final value u ^R (θ) of the Mth iteration; u ^R (θ)=F ^-1 (u ^R* (θ)) (1); Wherein u ^R (θ) in formula (1) represents the inverse Fourier transform result of the weighted frequency domain data u ^R* (θ) under the angle θ, and F ^-1 represents the inverse Fourier transform; The seventh step is to determine whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step; The eighth step is to use the final value u ^R (θ) of the Mth iteration as the initialization projection, let N=N+1, and enter the third step; The ninth step, specifically, takes the final value u ^R (θ) of the Mth iteration as the final recovered data result. 3. the recovery method of incomplete CT projection data under circular orbit according to claim 2, is characterized in that, The first step, obtaining the bad pixel position information is to convert the actually collected projection data into the corresponding chord diagram, that is: the x-axis of the chord diagram is the projected pixel position, and the y-axis is the projection angle position corresponding to the pixel position. ; Under the condition that a certain detector unit is damaged, a vertical line will appear on the corresponding pixel position in the chord diagram. According to this theory, the position information of the bad pixel point is obtained in the chord diagram, and the position of the bad pixel point is recorded. .

Images