Rapid photon counting CT projection data correction method

Abstract

The invention discloses a rapid photon counting CT projection data correction method. The rapid photon counting CT projection data correction method comprises the following steps that S1, a detector collects multi-energy-level projection data of a calibration part; s2, approximately averaging the multi-energy-level projection data into quasi-monoenergetic projection data; s3, looking up the table to obtain image domain data of the calibration piece; s4, converting the image domain data into projection domain data, and making up data loss caused by Radon transformation through an iteration mode; s5, acquiring projection data and projection data for table look-up; and S6, correcting projection data. In the calibration method provided by the invention, in order to make up for the limitations of high acquisition cost and large calculation amount of a single energy source in the existing method, a material attenuation coefficient, namely image domain data, under single energy is obtained by adopting a table look-up mode; image domain data obtained through table look-up are converted into projection domain data through Radon conversion, and data loss caused by Radon conversion is made up in an iteration mode.

Description

Technical Field

[0001] The invention belongs to the technical field of medical X-ray imaging, and in particular is a fast photon counting CT projection data correction method. Background Art

[0002] Computed tomography (CT) technology has become one of the most important clinical medical imaging methods due to its fast scanning speed, high spatial resolution and wide range of applications. With the widespread application of CT examinations, the limitations of traditional CT examinations such as high radiation dose and low contrast have become increasingly prominent. Photon counting spectral CT is based on the principle that the attenuation of X-rays in different substances is energy-dependent. It can provide specific tissue imaging and contrast agent concentration distribution, improve CT spatial resolution, enhance the contrast between lesions and healthy tissues, and reduce radiation dose. It is considered to be the next generation of CT technology.

[0003] However, the current photon counting detector technology is still imperfect. For example, due to limited signal processing capabilities, pulse accumulation will occur, resulting in a decrease in the number of collected photons and an increase in energy. Also, due to the inability to collect charges completely, charge sharing will occur, resulting in an increase in the number of photons and a decrease in energy. In addition, there are non-ideal factors such as x-ray escape and reabsorption, Compton scattering, and Rayleigh scattering, which cause the actual collected projection data to be inconsistent with the theoretical value, which will ultimately seriously affect the decomposition of the material and reduce the quality of the reconstructed image.

[0004] How to correct the projection data errors caused by non-ideal factors such as pulse accumulation, charge sharing, x-ray escape and reabsorption, Compton scattering, Rayleigh scattering, etc., to avoid affecting the subsequent material decomposition and image reconstruction is an urgent problem to be solved in the existing technology. Existing technical solutions often have certain limitations. For example, in the calibration scheme using a specific radiation source, the cost of the specific radiation source is very high. In addition, the more single the energy of the specific radiation source, the more accurate the calibration result, but in fact it is often difficult to obtain a radiation source with a single energy. In the calibration scheme using physical modeling, not only a deep physical background is required, and familiarity with various physical reactions in the acquisition process, but also the amount of calculation of physical modeling is very large. The classic Monte Carlo simulation modeling takes several days to run once. In the calibration scheme using a lookup table, the calibration process is limited to the image domain. In addition, in the process of image domain calibration, multiple reconstructions are often required, which is difficult to meet the real-time requirements in the actual scanning process.

[0005] Therefore, it is necessary to propose a fast and effective photon counting CT projection data correction method to reduce the calibration cost, reduce the calibration calculation amount, and meet the real-time requirements in the actual scanning process.

[0006] Non-ideal factors such as pulse accumulation, charge sharing, x-ray escape and reabsorption, Compton scattering, Rayleigh scattering, etc. can cause projection data errors. To address the above problems:

[0007] Patent CN 112415568A proposes a method of obtaining the particle energy measured by a specific detection channel by irradiating the detector with a specific radiation source, calculating a correction coefficient based on the obtained particle energy and the actual particle energy, and then performing energy correction based on the correction coefficient.

[0008] Patent CN 116011327A proposes a method of building a physical model, generating a data set based on the physical model, and using it to train a neural network. Then, a calibration experiment is performed to reduce the error between the neural network prediction value and the actual response value, and the physical parameters of the CT detector are determined.

[0009] CN 104090329A patent proposes a method of setting a fixed calibration part, storing the theoretical projection data of the fixed calibration part, and calculating the corresponding theoretical reconstructed image. Then, the actual projection data of the fixed calibration part is collected, and the actual reconstructed image is calculated. Finally, a mapping function is established to correct the actual reconstructed image to the theoretical value.

[0010] Patent CN 112415569A proposes a method of using a target detector and a reference detector to collect multiple samples of the same specific radiation source, and selecting at least one pair of completely deposited scattering events from the multiple samples, and then performing energy correction on the target detector based on the particle energy collected by the target detector, the particle energy collected by the reference detector and the actual particle energy.

[0011] Non-ideal factors such as pulse accumulation, charge sharing, x-ray escape and reabsorption, Compton scattering, and Rayleigh scattering may lead to errors in projection data, which will seriously affect the accuracy of subsequent material decomposition, reduce the quality of reconstructed images, and ultimately affect the physician's clinical diagnosis.

[0012] Existing calibration technologies can be divided into three categories: using specific radiation sources, physical modeling, and image domain calibration. The method using specific radiation sources has the limitations that the specific radiation sources are expensive and single-energy radiation sources are difficult to obtain. The physical modeling calibration method has the limitations that it requires familiarity with the physical process of the interaction between photons and detectors and has a huge amount of calculation. The image domain calibration method has the limitation that it is difficult to meet real-time requirements. Summary of the invention

[0013] In order to solve the defects and shortcomings in the above-mentioned prior art, the present invention provides a fast photon counting CT projection data correction method which can reduce calibration cost, reduce calibration calculation amount and meet the real-time requirements in the actual scanning process.

[0014] The technical solution adopted by the present invention to solve the technical problem is: a fast photon counting CT projection data correction method, comprising the following steps:

[0015] S1. The detector collects multi-level projection data of the calibration part;

[0016] S2. Approximately averaging the multi-energy level projection data into quasi-monoenergetic projection data;

[0017] S3. Obtain the single energy attenuation coefficient of the calibration piece by looking up the table in NI ST, that is, the image domain data of the calibration piece;

[0018] S4. Use Radon transform to convert the image domain data obtained by table lookup into projection domain data, and compensate for the data loss caused by Radon transform through iteration;

[0019] S5. Polynomial fitting of collected projection data and table lookup of projection data;

[0020] S6. Correct the projection data using the polynomial coefficients.

[0021] Preferably, step S1 specifically refers to: using a uniform material object of known material type as a calibration piece, and using a detector to collect multi-level projection data of the calibration piece; using a uniform material object of known material type is to facilitate subsequent table lookup to obtain the attenuation coefficient of the material.

[0022] Preferably, step S2 specifically means that: since the multi-level projection data collected by the photon detector is projection data under multiple energy window ranges, in order to facilitate comparison with the projection data under a single energy obtained by subsequent table lookup conversion, it is necessary to approximate the average of the projection data under the energy window range to be the projection data under a quasi-single energy, and the calculation method of the approximate average is:

[0023] Assuming that the energy window of a certain energy level is set to [i1,i2], then the quasi-monoenergetic E i =i1+(i2-i1)*0.5.

[0024] Preferably, step S3 specifically means: the NIST standard provides the attenuation coefficient of photons of different energies passing through a certain material, and the monoenergetic attenuation coefficient of a set of calibration parts is obtained by table lookup, and the attenuation coefficient obtained by table lookup is used as the ideal attenuation coefficient, that is, the ideal image domain data.

[0025] Preferably, the conversion between the projection domain data and the image domain data in step S4 mainly relies on Radon transform and inverse Radon transform; Radon transform and inverse Radon transform can be expressed as:

[0026] Let R be the Radon transform operator, then the Radon transform can be defined as:

[0027]

[0028] The inverse Radon transform can be obtained by Fourier integrating the variable S of R(S, α):

[0029]

[0030] make

[0031] Then perform a two-dimensional inverse Fourier transform to obtain the inverse Radon transform:

[0032]

[0033] From a mathematical analysis, the essence of Radon transform and inverse Radon transform is integral transform and Fourier transform, and both integral transform and Fourier transform are reversible, that is, in theory, ideal projection domain data can be obtained by ideal image domain data transformation. But in fact, taking the most commonly used implementation method of Radon inverse transform, the filtered back projection algorithm, as an example, since interpolation operations and Gibbs phenomenon will introduce errors in the transformation process, it is necessary to reduce the errors through iteration to compensate for the data loss caused by Radon transform.

[0034] Preferably, the image domain data obtained by table lookup is converted into projection domain data by using Radon transform, and the data loss caused by Radon transform is compensated by iteration; the specific steps are as follows: first, the monoenergetic attenuation coefficient of the known material, that is, the ideal image domain data, is obtained by table lookup, and the projection domain data is obtained by using Radon transform, and then the image domain data is obtained by using inverse Radon transform, and the error between the ideal image domain data obtained by table lookup and the image domain data obtained by the transformation is calculated, and the root mean square error RMSE indicator is used for the calculation error, and the specific calculation method is:

[0035]

[0036] Where N is the number of data in the ideal image domain data, x i is the ideal image domain data value, x i ′ is the image domain data value obtained by transformation.

[0037] Preferably, it is finally determined whether the error meets the requirement. In order not to affect the subsequent calculation of the CT value, the error needs to be controlled within one ten-thousandth of the maximum value of the ideal image domain data, that is, one ten-thousandth of the maximum value of the ideal image domain data is used as the threshold, and it is determined whether the calculated RMSE value is less than the threshold;

[0038] If it is greater than the threshold, it is necessary to update the projection domain data and then calculate the error, and determine whether it is less than the threshold. The specific update method is:

[0039] I′=I-radon(err)*k

[0040] Among them, I′ is the updated projection domain data, I is the projection domain data before the update, radon

[0041] represents Radon transform, which transforms image domain data back to projection domain data. err is the difference between the ideal image domain data and the image domain data obtained by the last updated projection domain data through Radon inverse transform, and k represents the transformation coefficient;

[0042] If it is less than the threshold, the transformed projection domain data or the last updated projection domain data is regarded as the ideal projection domain data.

[0043] Preferably, step S5 specifically refers to: using a polynomial fitting method to fit the collected quasi-monoenergetic projection data with the ideal projection data obtained by the above table lookup, and the specific fitting method is:

[0044]

[0045] Among them, a i is the polynomial coefficient, x is the projection domain data to be corrected, and y is the ideal projection domain data. The least squares method is used to solve the polynomial coefficient:

[0046]

[0047] Preferably, step S6 specifically refers to: using the polynomial coefficient a obtained in steps S1-S5 i And the collected projection data is used to obtain the corrected projection domain data.

[0048] In the calibration method proposed in the present invention, in order to make up for the limitations of the existing methods of high cost and large amount of calculation for obtaining a single energy source, a table lookup method is used to obtain the material attenuation coefficient under a single energy, that is, image domain data; then, the Radon transform is used to convert the image domain data obtained by the table lookup into projection domain data, and the data loss caused by the Radon transform is compensated by an iterative method; secondly, since the data obtained by the table lookup is single-energy data, and the actually collected data is multi-energy level data, the multi-energy level data is converted into quasi-monoenergetic data by approximate averaging; finally, a polynomial is used to fit the actually collected data and the data obtained by the table lookup, and the projection domain data is corrected according to the polynomial coefficients. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 It is a schematic diagram of the projection data correction process of the present invention;

[0050] Figure 2 It is a schematic diagram of the data flow of calculating the ideal projection domain of the present invention;

[0051] Figure 3 It is a relationship diagram between the number of iterations and the RMSE value of the present invention. DETAILED DESCRIPTION

[0052] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0053] Example 1

[0054] A fast photon counting CT projection data correction method, the basic process is as follows Figure 1 As shown, the method mainly includes the following steps:

[0055] The detector collects multi-level projection data of the calibration part: a uniform material object of known material type is used as the calibration part, and the detector is used to collect multi-level projection data of the calibration part. The uniform material object of known material type is used to facilitate the subsequent table lookup to obtain the attenuation coefficient of the material.

[0056] Approximately average the multi-energy level projection data into quasi-single-energy projection data: Since the multi-energy level projection data collected by the photon detector is projection data under multiple energy windows, in order to facilitate the comparison with the projection data under a single energy obtained by subsequent table lookup conversion, it is necessary to approximate the projection data under the energy window range to average the projection data under a quasi-single energy. The calculation method of the approximate average is:

[0057] Assuming that the energy window of a certain energy level is set to [i1,i2], then the quasi-monoenergetic E i =i1+(i2-i1)*0.5.

[0058] The monoenergetic attenuation coefficient of the calibration part is obtained by looking up the table in NI ST, that is, the image domain data of the calibration part: The NI ST standard provides the attenuation coefficient of photons of different energies passing through a certain material. Therefore, the monoenergetic attenuation coefficient of a set of calibration parts can be obtained by looking up the table, and the attenuation coefficient obtained by looking up the table is used as the ideal attenuation coefficient, that is, the ideal image domain data.

[0059] Convert the attenuation coefficient obtained by table lookup into projection domain data: The conversion between projection domain data and image domain data mainly relies on Radon transform and inverse Radon transform. Radon transform and inverse Radon transform can be expressed as:

[0060] Let R be the Radon transform operator, then the Radon transform can be defined as:

[0061]

[0062] The inverse Radon transform can be obtained by Fourier integrating the variable S of R(S, α):

[0063]

[0064] make

[0065] Then perform a two-dimensional inverse Fourier transform to obtain the inverse Radon transform:

[0066]

[0067] From a mathematical analysis, the essence of Radon transform and inverse Radon transform is integral transform and Fourier transform, and both integral transform and Fourier transform are reversible, that is, in theory, ideal projection domain data can be obtained by ideal image domain data transformation. But in fact, taking the most commonly used implementation method of Radon inverse transform, the filtered back projection algorithm, as an example, since interpolation operations and Gibbs phenomenon will introduce errors in the transformation process, it is necessary to reduce the errors through iteration to compensate for the data loss caused by Radon transform.

[0068] Based on this, Radon transform can be used to convert the image domain data obtained by table lookup into projection domain data, and the data loss caused by Radon transform can be compensated by iteration. The specific steps are as follows Figure 2 As shown: First, the single energy attenuation coefficient of the known material is obtained by looking up the table, that is, the ideal image domain data, and the projection domain data is obtained by Radon transform. Then, the image domain data is obtained by inverse Radon transform, and the error between the ideal image domain data obtained by looking up the table and the image domain data obtained by the transformation is calculated. The error is calculated using the root mean square error RMSE indicator. The specific calculation method is:

[0069]

[0070] Where N is the number of data in the ideal image domain data, x i is the ideal image domain data value, x i ′ is the image domain data value obtained by transformation.

[0071] Finally, it is determined whether the error meets the requirements. In order not to affect the subsequent calculation of the CT value, the error needs to be controlled within one ten-thousandth of the maximum value of the ideal image domain data. That is, one ten-thousandth of the maximum value of the ideal image domain data is used as the threshold to determine whether the calculated RMSE value is less than the threshold.

[0072] If it is greater than the threshold, it is necessary to update the projection domain data and then calculate the error, and determine whether it is less than the threshold. The specific update method is:

[0073] I′=I-radon(err)*k

[0074] Among them, I′ is the updated projection domain data, I is the projection domain data before the update, radon

[0075] represents Radon transform, which transforms image domain data back to projection domain data. err is the difference between the ideal image domain data and the image domain data obtained by the last updated projection domain data through inverse Radon transform, and k represents the transform coefficient.

[0076] If it is less than the threshold, the transformed projection domain data or the last updated projection domain data is regarded as the ideal projection domain data.

[0077] Fitting the collected projection data and the projection data obtained by table lookup: The collected quasi-monoenergetic projection data and the ideal projection data obtained by the table lookup are fitted by polynomial fitting. The specific fitting method is:

[0078]

[0079] Among them, a i is the polynomial coefficient, x is the projection domain data to be corrected, and y is the ideal projection domain data. The least squares method is used to solve the polynomial coefficient:

[0080]

[0081] Use polynomial coefficients to correct projection data: Use the polynomial coefficients a obtained in the above steps i And the collected projection data is used to obtain the corrected projection domain data.

[0082] Example 2

[0083] Assuming that a uniform Cu material is used as the calibration piece, the four energy level projection data of the calibration piece collected by the detector are as follows:

[0084] 10kev to 30kev energy range: [000.0000 000.0000 000.0000 000.0000] [113.2500 113.1638 112.9051 112.4740] [792.7500 792.7672 792.8190 792.9052] [906.0000 906.1380 906.5519 907.2416] [792.7500 792.7672 792.8190 792.9052] [113.2500 113.1638 112.9051 112.4740] [000.0000 000.0000 000.0000 000.0000]

[0092] 30kev to 50kev energy range: [000.0000 000.0000 000.0000 000.0000] [016.1250 016.1127 016.0759 016.0145] [112.8750 112.8775 112.8848 112.8971] [129.0000 129.0196 129.0786 129.1768] [112.1250 112.8775 112.8848 112.8971] [016.1250 016.1127 016.0759 016.0145] [000.0000 000.0000 000.0000 000.0000]

[0100] 50kev to 70kev energy range: [000.0000 000.0000 000.0000 000.0000] [005.2500 005.2460 005.2340 005.2140] [036.7500 036.7508 036.7532 036.7572] [042.0000 042.0064 042.0256 042.0576] [036.7500 036.7508 036.7532 036.7572] [005.2500 005.2460 005.2340 005.2140] [000.0000 000.0000 000.0000 000.0000]

[0108] 70kev to 90kev energy range: [000.0000 000.0000 000.0000 000.0000] [002.6250 002.6230 002.6170 002.6070] [018.3750 018.3754 018.3766 018.3786] [021.0000 021.0032 021.0128 021.0288] [018.3750 018.3754 018.3766 018.3786] [002.6250 002.6230 002.6170 002.6070] [000.0000 000.0000 000.0000 000.0000]

[0116] Convert the multi-energy projection data into quasi-mono-energy projection data: According to the calculation method of quasi-mono-energy, then: E1 = 10 + (30-10) * 0.5 = 20, E2 = 30 + (50-30) * 0.5 = 40, E3 =

[0117] 50+(70-50)*0.5=60, E4=70+(90-70)*0.5=80.

[0118] The monoenergetic attenuation coefficient of a specific material is obtained by looking up the table in NIST: the mass attenuation coefficients of the Cu element at 20kev, 40kev, 60kev, and 80kev are 33.79, 4.862, 1.593, and 0.763, respectively. The attenuation coefficient obtained by looking up the table is the mass attenuation coefficient of the material, while the CT image domain data is the linear attenuation coefficient. The conversion relationship between the mass attenuation coefficient and the linear attenuation coefficient is as follows:

[0119] μ x =μ*ρ x

[0120] Among them, μ x is the linear attenuation coefficient, μ is the mass attenuation coefficient, ρ x is the density of the material. It is known that the density of Cu element is 8.96. Substituting it into the above formula, we can obtain the linear attenuation coefficients of Cu element at 20kev, 40kev, 60kev and 80kev energies are 302.7584, 43.56352, 14.27328 and 6.83648 respectively.

[0121] Convert the attenuation coefficient obtained from the table into projection domain data:

[0122] Taking 20kev energy as an example, since the calibration piece used is a uniform Cu element, and the physical meaning of the image domain data value is the attenuation coefficient of the material passing through the point, therefore, the ideal image domain data value should be the attenuation coefficient of the Cu element at 20kev energy, that is, 302.7584.

[0123] Radon transform to projection domain data is: [000.0000 000.0000 000.0000 000.0000] [113.5344 113.4479 113.1886 112.7564] [794.7408 794.7581 794.8100 794.8964] [908.2752 908.4135 908.8285 909.5200] [794.7408 794.7581 794.8100 79.8964] [113.5344 113.4479 113.1886 112.7564] [000.0000 000.0000 000.0000 000.0000]

[0131] The inverse Radon transform to image domain data is: [292.5996 210.0267 282.4276] [294.8746 207.7517 294.8746] [282.4276 210.0267 292.5996]

[0135] The error between the ideal image domain data obtained by table lookup and the image domain data obtained by transformation is:

[0136] [-10.1588 -92.7317 -20.3308]

[0137] [-07.8838 -95.0067 -07.8838]

[0138] [-20.3308 -92.7317 -10.1588]

[0139] The RMSE value is 55.1584, and the threshold is set to one ten-thousandth of the maximum value of the ideal image domain data, that is, 302.7584*0.0001=0.0302. If the RMSE value is greater than the threshold, the projection domain data is updated.

[0140] The error matrix is ​​converted into projection domain data: [0000.0000 0000.0000 0000.0000 0000.0000]

[0142] [-004.7967 -004.8818 -004.9596 -005.0300]

[0143] [-063.8389 -063.8350 -063.8234 -063.8040]

[0144] [-219.9460 -219.7836 -219.6513 -219.5492]

[0145] [-063.8389 -063.8350 -063.8234 -063.8040]

[0146] [-004.7967 -004.8818 -004.9596 -005.0300] [0000.0000 0000.0000 0000.0000 0000.0000]

[0148] The updated projection domain data is: [0000.0000 0000.0000 0000.0000 0000.0000] [0118.3311 0118.3297 0118.1482 0117.7864] [0858.5797 0858.5931 0858.6333 0858.7004] [1128.2213 1128.1974 1128.4800 1129.0694] [0858.5797 0858.5931 0858.6333 0858.7004] [0118.3311 0118.3297 0118.1482 0117.7864] [0000.0000 0000.0000 0000.0000 0000.0000]

[0156] After the update, the projection domain data is transformed back to the image domain, and the RMSE value between it and the ideal image domain data is 34.2982.

[0157] The RMSE value after 100 iterations is as follows Figure 3 As shown. Among them, after 87 iterations, the RMSE value is less than the threshold value of 0.0302, meeting the error requirement. At this time, the ideal projection domain data is: [0000.0000 0000.0000 0000.0000 0000.0000] [0126.3549 0131.9814 0137.4136 0142.6501] [0905.2371 0905.2532 0905.3014 0905.3818] [1135.3562 1124.0700 1113.1095 1102.4763] [0905.2371 0905.2532 0905.3014 0905.3818] [0126.3549 0131.9814 0137.4136 0142.6501] [0000.0000 0000.0000 0000.0000 0000.0000]

[0165] Polynomial fitting of the collected projection data and the projection data from the table: the least square method is used to fit the polynomial coefficients. For example, a 10th-order polynomial fitting is used, and the polynomial coefficients are:

[0166] [-0.0000 0.0000-0.0000 0.0000-0.0000 0.0000-0.0000 0.0023-0.17245.1146 0.0000]

[0167] Use polynomial coefficients to correct projection data: The projection domain data after correction using polynomial coefficients is: [0000.0000 0000.0000 0000.0000 0000.0000] [0126.2214 0132.1791 0137.3345 0142.6633] [0905.3844 0905.2463 0905.0566 0905.4964] [1135.6350 1126.6620 1112.0960 1102.6530] [0905.3844 0905.2463 0905.0566 0905.4964] [0126.2214 0132.1791 0137.3345 0142.6633]

[0174] [0000.0000 0000.0000 0000.0000 0000.0000].

[0175] Example 3

[0176] To prove the effectiveness of the present invention, the projection domain data before and after correction are converted into image domain data, and the difference with the ideal image domain data is calculated. Peak Signal to Noise Ratio (PSNR) and Mean Absolute Error (MAE) indicators are used for comparison. The average value is taken after 10 tests. The calculation results are shown in Table 1. The results show that by using the correction method proposed by the present invention, the corrected projection data is significantly improved compared with the projection data before correction. In addition, when the data size is 512*512, the method of the present invention only takes 15.0279 seconds to calculate the polynomial coefficients. During the real-time inspection process, it only takes 0.8145 seconds to correct the projection data using the polynomial coefficients, which meets the real-time requirements.

[0177] Table 1 Differences between the ideal image domain data before and after correction

[0178] Before calibration After correction PSNR 13.2342 34.0483 MAE 40.3498 3.9493

[0179] In summary, the effectiveness of a fast photon technology CT projection data correction method provided by the present invention is proved. Although the specific embodiments of the present invention have been shown and described, it is understood by those skilled in the art that various changes, modifications, substitutions and variations can be made to these specific embodiments without departing from the principles and spirits, and the scope of the present invention is defined by the attached claims and their equivalents.

Claims

1. A fast photon counting CT projection data correction method, characterized in that: The following steps are involved: S1. The detector collects multi-level projection data of the calibration part; S2. Approximately averaging the multi-energy level projection data into quasi-monoenergetic projection data; S3. Obtain the monoenergetic attenuation coefficient of the calibration piece by looking up the table in NIST, that is, the image domain data of the calibration piece; S4. Use Radon transform to convert the image domain data obtained by table lookup into projection domain data, and compensate for the data loss caused by Radon transform through iteration; S5. Polynomial fitting of the collected projection data and the projection data obtained by table lookup; S6. Correct the projection data using the polynomial coefficients.

2. A fast photon counting CT projection data correction method according to claim 1, characterized in that: The step S1 specifically refers to: using a uniform material object of known material type as a calibration piece, and using a detector to collect multi-level projection data of the calibration piece; using a uniform material object of known material type is to facilitate subsequent table lookup to obtain the attenuation coefficient of the material.

3. A fast photon counting CT projection data correction method according to claim 1, characterized in that: The step S2 specifically means that since the multi-level projection data collected by the photon detector is projection data under multiple energy window ranges, in order to facilitate comparison with the projection data under a single energy obtained by subsequent table lookup conversion, it is necessary to approximate the average of the projection data under the energy window range to be the projection data under a quasi-single energy. The calculation method of the approximate average is: Assuming that the energy window of a certain energy level is set to [i1,i2], then the quasi-monoenergetic E i =i1+(i2-i1)*0.

5.

4. The method for correcting fast photon counting CT projection data according to claim 1, characterized in that: The step S3 specifically means that the NIST standard provides the attenuation coefficient of photons of different energies passing through a certain material, and the monoenergetic attenuation coefficient of a set of calibration parts is obtained by table lookup, and the attenuation coefficient obtained by table lookup is used as the ideal attenuation coefficient, that is, the ideal image domain data.

5. The method for correcting fast photon counting CT projection data according to claim 1, characterized in that: The conversion between the projection domain data and the image domain data in step S4 mainly relies on Radon transform and inverse Radon transform; Radon transform and inverse Radon transform can be expressed as: Let R be the Radon transform operator, then the Radon transform can be defined as: The inverse Radon transform can be obtained by Fourier integrating the variable S of R(S, α): make Then perform a two-dimensional inverse Fourier transform to obtain the inverse Radon transform: From a mathematical analysis, the essence of Radon transform and inverse Radon transform is integral transform and Fourier transform, and both integral transform and Fourier transform are reversible, that is, in theory, ideal projection domain data can be obtained by ideal image domain data transformation; but in fact, taking the most commonly used implementation method of inverse Radon transform, the filtered back projection algorithm, as an example, since interpolation operations and Gibbs phenomenon will introduce errors in the transformation process, it is necessary to reduce the errors through iteration to compensate for the data loss caused by Radon transform.

6. A fast photon counting CT projection data correction method according to claim 5, characterized in that: The Radon transform is used to convert the image domain data obtained by table lookup into projection domain data, and the data loss caused by Radon transform is compensated by iteration. The specific steps are as follows: first, the monoenergetic attenuation coefficient of the known material, that is, the ideal image domain data, is obtained by table lookup, and the projection domain data is obtained by Radon transform. Then, the image domain data is obtained by inverse Radon transform, and the error between the ideal image domain data obtained by table lookup and the image domain data obtained by transformation is calculated. The root mean square error RMSE indicator is used for the calculation error. The specific calculation method is: Where N is the number of data in the ideal image domain data, x i is the ideal image domain data value, x i ′ is the image domain data value obtained by transformation.

7. A fast photon counting CT projection data correction method according to claim 6, characterized in that: Finally, it is determined whether the error meets the requirements. In order not to affect the subsequent calculation of CT values, the error needs to be controlled within one ten-thousandth of the maximum value of the ideal image domain data. That is, one ten-thousandth of the maximum value of the ideal image domain data is used as the threshold to determine whether the calculated RMSE value is less than the threshold. If it is greater than the threshold, it is necessary to update the projection domain data and then calculate the error, and determine whether it is less than the threshold. The specific update method is: I′=I-radon(err)*k Where I′ is the updated projection domain data, I is the projection domain data before the update, radon represents Radon transform, which transforms the image domain data back to the projection domain data; err is the difference between the ideal image domain data and the image domain data obtained by the last updated projection domain data through Radon inverse transform, and k represents the transformation coefficient; If it is less than the threshold, the transformed projection domain data or the last updated projection domain data is regarded as the ideal projection domain data.

8. The method for correcting fast photon counting CT projection data according to claim 1, characterized in that: Step S5 specifically refers to: using a polynomial fitting method to fit the collected quasi-monoenergetic projection data with the ideal projection data obtained by the above table lookup, and the specific fitting method is: Among them, a i is the polynomial coefficient, x is the projection domain data to be corrected, and y is the ideal projection domain data; the least squares method is used to solve the polynomial coefficient:

9. The method for correcting fast photon counting CT projection data according to claim 1, characterized in that: Step S6 specifically refers to: using the polynomial coefficients ai solved in steps S1-S5 and the collected projection data to obtain corrected projection domain data.

Images