A dual-energy spectrum material decomposition calibration phantom and method

Abstract

The invention relates to the field of CT technology, and discloses a dual-energy spectrum material decomposition calibration phantom and method, comprising a 16 cm long organic glass cylinder, and two Teflon cylinders with different diameters are inlaid on the organic glass cylinder, and the diameters of the two Teflon cylinders are respectively one-third and one-quarter of the diameter of the organic glass cylinder. The dual-energy spectrum material decomposition calibration phantom and method are calibrated on the projection, do not rely on the image reconstruction process, can make fuller use of the energy spectrum and the original projection data information, and obtain more accurate images and material decomposition results.

Description

Technical Field

[0001] The present invention relates to the field of CT technology, and in particular to a dual-energy spectrum material decomposition calibration phantom and method. Background Art

[0002] General CT has only one set of energy spectra. After scanning, the grayscale value of the image corresponds to the linear attenuation coefficient at a certain energy, or the CT value converted from it. Such CT equipment has no material resolution capability. Energy spectrum CT refers to a device that obtains scanning data under multiple energy spectra at the same time, and reconstructs the corresponding multi-energy spectra of these data to obtain images with material resolution.

[0003] Spectral calibration of spectral CT is crucial to obtain accurate quantitative scan images and physical resolution. A common approach is to calibrate and reconstruct spectral CT as multiple single-energy scans to obtain images at two energies, and then perform material decomposition on the two sets of images.

[0004] The common dual-energy decomposition method is to first obtain two sets of single-energy images by traditional methods, and then perform decomposition based on the images. This method is relatively simple and intuitive, but it is restricted by image quality. If there are artifacts on the image or the CT value is inaccurate, the material decomposition result will be greatly reduced. In addition, errors are easily generated during the calculation process, affecting the calculation results. For this reason, the present invention provides a dual-energy spectrum material decomposition calibration phantom and method. Summary of the invention

[0005] In view of the above-mentioned deficiencies in the prior art, the present invention provides a dual-energy spectrum material decomposition calibration phantom and method, which can reduce errors, improve calculation accuracy, and obtain more accurate images and material decomposition results.

[0006] The present invention provides the following technical solution: a dual-energy spectrum material decomposition calibration phantom, comprising a 16 cm long organic glass cylinder, and two Teflon cylinders with different diameters are inlaid on the organic glass cylinder, and the diameters of the two Teflon cylinders are respectively one third and one quarter of the diameter of the organic glass cylinder.

[0007] A dual-energy spectrum material decomposition calibration method uses a dual-energy spectrum material decomposition calibration phantom placed at different positions, performs several scans, and then directly calculates two sets of spectra of the X-ray source on the obtained projection data. The specific operation is as follows:

[0008] S1. The emission spectrum of an X-ray tube at one kVp is expressed as the superposition of the characteristic peak and the Gaussian peak. The formula is:

[0009] ;

[0010] Among them, tri() represents a triangular peak, E0 is the characteristic peak position of the X-ray tube target material, W0 is the peak width, A0 is the peak amplitude, and the second term represents the superposition of n Gaussian peaks, {A n , E n and W n}n>=1 is the amplitude, peak position and peak width of the Gaussian peak;

[0011] S2, according to the data P0 obtained by scanning the phantom, extract two tracks of the edge of the phantom from the sine graph, and then obtain the coordinates of the center of the phantom;

[0012] S3, the projection data P0 obtained by actual scanning and the theoretical projection data P are obtained by the following formula {A n , E n and W n}n>=1, thus obtaining the spectrum of X-rays, and the formula is:

[0013] {A n , E n , W n}=argmin(Projection(A n , E n , W n ; d1, d2)-P0) 2 ;

[0014] S4. Repeat this process for the data scanned at different energies to obtain two sets of spectra at different energies, calculate the projection data P, and use the formula P=Projection(An, En and Wn; d1, d2) to obtain two projection data at different energies, which are represented by P1 and P2 respectively;

[0015] S5, sorting the theoretical projection data at different angles according to the length of the organic glass and Teflon passed through on the ray path, to obtain a data pair of {path length d1 of the ray passing through the organic glass, path length d2 of the ray passing through the Teflon}=>{projection value P1, projection value P2};

[0016] S6. Use two-dimensional interpolation to upsample the above data pairs to obtain projection data corresponding to equally spaced combinations of plexiglass lengths and Teflon lengths within a certain range. This set of data is recorded as R: {d1, d2}=>{P1, P2};

[0017] S7. Perform corresponding air correction on the actual collected model data. After the detector is geometrically corrected, the obtained data is recorded as {P1', P2'}. Based on {P1', P2'}, a table lookup operation is performed in R to obtain the ideal material length {d1', d2'} corresponding to {P1', P2'};

[0018] S8. Use a nonlinear model F to fit the mapping relationship from {P1, P2} to {d1, d2}. Use the nonlinear least squares method to fit the data scanned in S4 to obtain the ten coefficients required by model F. The specific form of the nonlinear model is as follows:

[0019] ;

[0020] S9. In the process, the range of {d1, d2} should be ≥ 0. The system scanning range is defined as an array of equally spaced distributions between 0 and 200 mm. d1 uses 21 points with an equal interval of 10 mm between 0 and 200 mm, and d2 uses 21 points with an equal interval of 3 mm between 0 and 60 mm for sampling for subsequent calculations.

[0021] S10, supplement the projection data of d2 between 0 and 60 mm corresponding to d1 between 0 and -50 mm, and supplement the projection data of d1 between 0 and 200 mm corresponding to d2 between 0 and -15 mm. When d1 or d2 is a negative number, it is treated as 0 and the projection data is supplemented to make the corresponding P1 and P2 results positive;

[0022] S11. After the coefficients are obtained, they are saved as correction vectors. When actual scanning data arrives, a set of length values ​​can be directly obtained from a set of projection values ​​using model F.

[0023] S12. Based on the length values ​​at all angles, the images of the two materials, plexiglass and Teflon, can be reconstructed. After reconstructing d1 and d2 respectively to obtain the image A(x, y) of plexiglass and the image B(x, y) of Teflon, the virtual monoenergetic image can be used to display images at different energies according to actual needs. The formula for any given energy E is:

[0024] ;

[0025] According to the definition of CT value, the attenuation coefficient image is further converted into CT value through the attenuation of water under this energy for diagnosis. When E is small, the obtained image has higher soft tissue contrast; when E is large, the metal artifacts on the image will be reduced.

[0026] Preferably, when the scanning light source position, the detector position, and the rotation center are known in S2, the path length {d1, d2} of each ray passing through the phantom can be further calculated by using the obtained phantom position. n , E n and W n}n>=1, the projection data P on each detector pixel can be obtained theoretically, and the calculation formula is P=Projection(A n, E n and W n ; d1, d2), where Projection represents the orthographic projection based on spectrum and geometry, and the projection data P is calculated when each ray passes through the plexiglass and Teflon.

[0027] Compared with the prior art, the present invention has the following beneficial effects:

[0028] The dual-energy spectrum material decomposition calibration model and method are calibrated on the projection, do not rely on the image reconstruction process, can make fuller use of the energy spectrum and the original projection data information, and obtain more accurate images and material decomposition results; in addition, the d value is prone to negative values ​​during the calculation process, but d represents the physical length and there is no negative number, so the data set is expanded to ensure that d is a positive number, thereby reducing experimental errors and improving experimental accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 This is a schematic diagram of the calibration phantom structure of the present invention;

[0030] Figure 2 This is the shape diagram of the emission spectrum of the bulb tube of the present invention:

[0031] Figure 3 It is the theoretical projection data diagram under two sets of spectra of the present invention;

[0032] Figure 4 This is a data diagram of d1 and d2 of the present invention;

[0033] Figure 5 is a truncation effect diagram of the data set of the present invention;

[0034] Figure 6 It is the projection diagram of the data after expansion of the present invention;

[0035] Figure 7 This is a calculation diagram of the positive values ​​of d1 and d2 of the present invention. DETAILED DESCRIPTION

[0036] In order to make the objectives, technical solutions and advantages of the embodiments of the present disclosure clearer, the technical solutions of the embodiments of the present disclosure will be clearly and completely described below in conjunction with the drawings of the embodiments of the present disclosure. In order to keep the following description of the embodiments of the present disclosure clear and concise, the present disclosure omits detailed descriptions of known functions and known components to avoid unnecessary confusion of the concepts of the present invention.

[0037] See also Figure 1 A dual-energy spectrum material decomposition calibration phantom includes a 16-cm-long organic glass cylinder, and two Teflon cylinders with different diameters are inlaid on the organic glass cylinder, and the diameters of the two Teflon cylinders are respectively one-third and one-quarter of the diameter of the organic glass cylinder.

[0038] A dual-energy spectrum material decomposition calibration method uses a dual-energy spectrum material decomposition calibration phantom placed at different positions, performs several scans, and then directly calculates two sets of spectra of the X-ray source on the obtained projection data. The specific operation is as follows:

[0039] S1. The emission spectrum of an X-ray tube at one kVp is represented as the superposition of a characteristic peak and a Gaussian peak, such as Figure 2 , the formula is:

[0040] ;

[0041] Among them, tri() represents a triangular peak, E0 is the characteristic peak position of the X-ray tube target material, W0 is the peak width, A0 is the peak amplitude, the second term represents the superposition of n Gaussian peaks, {An, En and Wn}n>=1 is the amplitude, peak position and peak width of the Gaussian peak, for the common tungsten target, It is taken as a fixed value of 69.525 keV. According to the experimental test results, n can be an integer between 3 and 6.

[0042] S2. According to the data P0 obtained by scanning the phantom, two trajectories of the edge of the phantom are extracted from the sine graph, and the coordinates of the center of the phantom are obtained accordingly.

[0043] When the scanning light source position, detector position, and rotation center are known, the path length {d1, d2} of each ray passing through the phantom can be further calculated by using the obtained phantom position. n , E n and W n}n>=1, the projection data P on each detector pixel can be obtained theoretically, and the calculation formula is P=Projection(A n , E n and W n ; d1, d2), where Projection represents the orthographic projection based on spectrum and geometry. When each ray passes through plexiglass and Teflon, the projection data P is obtained by the formula.

[0044] S3, the projection data P0 obtained by actual scanning and the theoretical projection data P are obtained by the following formula {A n , E n and W n}n>=1, thus obtaining the spectrum of X-rays, and the formula is:

[0045] {A n , E n , W n}=argmin(Projection(An , E n , W n ; d1, d2)-P0) 2 .

[0046] S4. Repeat this process for the data scanned at different energies to obtain two sets of spectra at different energies, calculate the projection data P, and use the formula P=Projection(An, En and Wn; d1, d2) to obtain two projection data at different energies, which are represented by P1 and P2 respectively.

[0047] S5. Sort the theoretical projection data at different angles according to the length of the plexiglass and Teflon passed through on the ray path to obtain a data pair of {path length d1 of the ray passing through the plexiglass, path length d2 of the ray passing through Teflon}=>{projection value P1, projection value P2}.

[0048] S6. Use two-dimensional interpolation to upsample the above data pairs to obtain projection data corresponding to equally spaced combinations of plexiglass lengths and Teflon lengths within a certain range. This set of data is recorded as R: {d1, d2}=>{P1, P2}.

[0049] S7. Perform corresponding air correction on the model data actually collected. After the detector geometric correction, the obtained data is recorded as {P1', P2'}. Based on {P1', P2'}, a table lookup operation is performed in R to obtain the ideal material length {d1', d2'} corresponding to {P1', P2'}.

[0050] S8. Use a nonlinear model F to fit the mapping relationship from {P1, P2} to {d1, d2}. Use the nonlinear least squares method to fit the data scanned in S4 to obtain the ten coefficients required by model F. The specific form of the nonlinear model is as follows:

[0051] .

[0052] S9. In the process, the range of {d1, d2} should be ≥ 0. The system scanning range is defined as an array of equally spaced distributions between 0 and 200 mm. {P1, P2} are theoretical projection data corresponding to different lengths of {d1, d2} under two sets of spectra, and the digital range is 0 to 1. The corresponding relationship between {d1, d2} and P under one set of spectra is as follows: Figure 3 , d1 uses 21 points with an equal interval of 10mm between 0 and 200mm, and d2 uses 21 points with an equal interval of 3mm between 0 and 60mm for sampling for subsequent calculations. The calculation data is as follows Figure 4 .

[0053] If such a data set is used directly for curve fitting, the truncation effect occurs at the position where d1 and d2 are 0, causing the fitting curve to fluctuate greatly around 0, and it is possible to obtain a result where d1 and d2 are negative numbers, such as Figure 5 .

[0054] S10, supplement the projection data between 0 and 60 mm corresponding to d1 being 0 to -50 mm, and the projection data between 0 and 200 mm corresponding to d2 being 0 to -15 mm. When d1 or d2 is a negative number, it is treated as 0 to complete the data, such as Figure 6 , so that the corresponding results of P1 and P2 are positive, such as Figure 7 , and also when there is system noise in p1 and p2, the fitting result can still be guaranteed to be positive.

[0055] S11. After the coefficients are obtained, they are saved as correction vectors. When actual scanning data arrives, a set of length values ​​can be directly obtained from a set of projection values ​​using model F.

[0056] S10. Based on the length values ​​at all angles, the images of the two materials, plexiglass and Teflon, can be reconstructed. After reconstructing d1 and d2 respectively to obtain the image A(x, y) of plexiglass and the image B(x, y) of Teflon, the virtual monoenergetic image can be used to display images at different energies according to actual needs. The formula for any given energy E is:

[0057] ;

[0058] According to the definition of CT value, the attenuation coefficient image is further converted into CT value through the attenuation of water under this energy for diagnosis. When E is small, the obtained image has higher soft tissue contrast; when E is large, the metal artifacts on the image will be reduced.

[0059] The above embodiments are only exemplary embodiments of the present invention and are not intended to limit the present invention. The protection scope of the present invention is defined by the claims. Those skilled in the art may make various modifications or equivalent substitutions to the present invention within the essence and protection scope of the present invention, and such modifications or equivalent substitutions shall also be deemed to fall within the protection scope of the present invention.

Claims

1. A dual-energy spectrum material decomposition calibration method, characterized in that: The calibration model is a 16 cm long organic glass cylinder, and two Teflon cylinders with different diameters are inlaid on the organic glass cylinder, and the diameters of the two Teflon cylinders are respectively one third and one quarter of the diameter of the organic glass cylinder; The calibration phantom is placed in different positions, scanned several times, and then two sets of spectra of the X-ray source are calculated directly on the obtained projection data. The specific operation is as follows: S1. The emission spectrum of an X-ray tube at one kVp is expressed as the superposition of the characteristic peak and the Gaussian peak. The formula is: ； Among them, tri() represents a triangular peak, E0 is the characteristic peak position of the X-ray tube target material, W0 is the peak width, A0 is the peak amplitude, and the second term represents the superposition of n Gaussian peaks, n>=1, A n , E n and W n is the amplitude, peak position and peak width of the Gaussian peak; S2, according to the data P0 obtained by scanning the phantom, two trajectories of the edge of the phantom are extracted from the sine graph, and the coordinates of the center of the phantom are calculated accordingly; S3, the projection data P0 obtained by actual scanning and the theoretical projection data P are obtained by the following formula {A n , E n , W n } to obtain the spectrum of X-rays, and the formula is: {AT n ,IT n ,W n }=argmin(Projection(A n ,IT n ,W n ;d1,d2)-P0) 2 ; S4. Repeat this process for the data scanned at different energies to obtain two sets of spectra at different energies, calculate the projection data P, and calculate the projection data at two different energies by the formula P=Projection(An, En, Wn; d1, d2), which are represented by P1 and P2 respectively; S5, sorting the theoretical projection data at different angles according to the length of the organic glass and Teflon passed through on the ray path, to obtain a data pair of {path length d1 of the ray passing through the organic glass, path length d2 of the ray passing through the Teflon}=>{projection value P1, projection value P2}; S6. Use two-dimensional interpolation to upsample the above data pairs to obtain projection data corresponding to equally spaced combinations of plexiglass lengths and Teflon lengths within a certain range. This set of data is recorded as R: {d1, d2}=>{P1, P2}; S7. Perform corresponding air correction on the actual collected model data. After the detector is geometrically corrected, the obtained data is recorded as {P1', P2'}. Based on {P1', P2'}, a table lookup operation is performed in R to obtain the ideal material length {d1', d2'} corresponding to {P1', P2'}; S8. Use a nonlinear model F to fit the mapping relationship from {P1, P2} to {d1, d2}. Use the nonlinear least squares method to fit the data scanned in S4 to obtain the ten coefficients required by model F. The specific form of the nonlinear model is as follows: ； S9, during the process, the range of {d1, d2} is greater than or equal to 0, the system scanning range is defined as an array of equally spaced distributions between 0 and 200 mm, d1 selects 21 points with an equal interval of 10 mm between 0 and 200 mm, and d2 selects 21 points with an equal interval of 3 mm between 0 and 60 mm for sampling for subsequent calculations; S10, supplement the projection data of d2 between 0 and 60 mm corresponding to d1 between 0 and -50 mm, and supplement the projection data of d1 between 0 and 200 mm corresponding to d2 between 0 and -15 mm, and when d1 or d2 is negative, treat it as 0 to complete the data, so that the corresponding results of P1 and P2 are positive values; S11, after the coefficients are obtained, they are saved as correction vectors. When the actual scanning data arrives, a set of length values ​​is directly obtained from a set of projection values ​​using model F; S12. Reconstruct the images of the two materials, plexiglass and Teflon, based on the length values ​​at all angles. After reconstructing d1 and d2 respectively to obtain the image A(x, y) of plexiglass and the image B(x, y) of Teflon, use the virtual monoenergetic image to display the images at different energies according to actual needs. The formula for any given energy E is: ； According to the definition of CT value, the attenuation coefficient image is further converted into CT value through the attenuation of water under this energy for diagnosis. When E is small, the obtained image has higher soft tissue contrast; when E is large, the metal artifacts on the image will be reduced.

2. A dual-energy spectrum material decomposition calibration method according to claim 1, characterized in that: When the scanning light source position, detector position, and rotation center are known in S2, the path length {d1, d2} of each ray passing through the phantom is further calculated by using the obtained phantom position. n , E n , W n }, the projection data P on each detector pixel is calculated theoretically, and the calculation formula is P=Projection(A n , E n , W n ; d1, d2), where Projection represents the orthographic projection based on spectrum and geometry, and the projection data P is obtained when each ray passes through plexiglass and Teflon.

Images