Method and computer program for generating channel spectrum response difference correction table and channel spacing correction table

Abstract

The present invention discloses a method for generating a channel spectrum response difference correction table and a channel spacing correction table, a computer program, a computer readable storage medium and a computer tomography device. The method for generating a channel spectrum response difference correction table and a channel spacing correction table includes using parallel projection P paral (k, r) determines the channel position corresponding to the center of the phantom, and calculates the simulated parallel projection and the ideal projection P accordingly sim (k, r), and finally according to the projection difference diff(k, l, r) and the simulated ideal projection P sim (k, r), calculate the channel spectrum response difference correction table t CCR (k, l) and channel spacing correction table t spacing The method, computer program, computer-readable storage medium, and computer tomography device of the present invention for generating a channel spectrum response difference correction table and a channel spacing correction table are not affected by residual beam hardening and can effectively generate ideal projections.

Description

[0001] Nouns and pronouns referring to persons in this patent application are not limited to a specific gender. Technical Field

[0002] The present invention relates to computed tomography, and in particular to channel spectral response difference correction and channel spacing correction. Background Art

[0003] The X-ray beam generated by an X-ray tube covers a wide frequency spectrum that varies with the increasing attenuation of the object. The varying spectral responses of the detector channels can cause signal offsets, resulting in ring artifacts in the image. Furthermore, the inevitable non-equidistant spacing of detector cells due to mechanical errors can cause ring and stripe artifacts in the image. To prevent these artifacts, channel spectral response difference correction (CCR) and channel spacing correction (spacing) are required.

[0004] The key to generating the channel spectral response difference correction table and the channel spacing correction table is to generate the ideal projection. The currently used fitting method has achieved good results. However, due to the inevitable residual beam hardening, the fitting method causes the fitting results to deviate from the actual ideal projection, especially when some projections are missing. Summary of the Invention

[0005] In view of this, the present invention proposes a method, a computer program, a computer-readable storage medium, and a computed tomography device for generating a channel spectrum response difference correction table and a channel spacing correction table.

[0006] According to a first aspect of the present invention, a method for generating a channel spectrum response difference correction table and a channel spacing correction table is provided, comprising: step S102, obtaining an original projection P(k, l, r) of a water phantom, wherein the water phantom can cover all channels of the detector, where k represents the channel number, l represents the sorting number, and r represents the reading number; step S104, averaging the original projections P(k, l, r) in the row direction to obtain an average projection p mean (k, r); Step S106, for each channel, interpolate the average projection P in the projection angle direction mean (k, r), to obtain the parallel projection P paral (k, r); Step S108, according to the parallel projection P paral (k, r), determine the channel position corresponding to the center of the phantom; step S110, calculate the simulated parallel projection according to the channel position corresponding to the center of the phantom Step S112: for each channel, interpolate and simulate parallel projection in the projection angle direction To obtain the simulated ideal projection P sim (k, r), and calculate the projection difference diff(k, l, r) between it and the original projection, diff(k, l, r) = P(k, l, r) - P sim (k, r); Step S114, according to the projection difference diff(k, l, r) and the simulated ideal projection P sim (k, r), calculate the channel spectrum response difference correction table t CCR (k, l) and channel spacing correction table t spacing (k, l).

[0007] In one embodiment, the step S108 includes: step S116, performing the parallel projection P paral (k, r) is convolution filtered; step S118, determining the channel positions E1(r) and E2(r) corresponding to the edges of the phantom in the filtered parallel projection; step S120, calculating the channel position C(r) corresponding to the center of the phantom based on the channel positions E1(r) and E2(r) corresponding to the edges of the phantom, C(r) = (E1(r) + E2(r)) / 2.

[0008] In one embodiment, step S116 includes performing convolution filtering according to the following formula:

[0009] FP=iFFT(FFT(p paral )×H(ω))

[0010] Among them, FFT() is Fourier transform, H(ω) is the convolution filter function, and iFFT() is the inverse Fourier transform.

[0011] In one embodiment, H(ω) is the product of a ramp function and a windowing function.

[0012] In one embodiment, the parallel projection P paral The missing data of (k, r) results in the loss of one of the channel positions E1(r) and E2(r) corresponding to the edge of the phantom. Assuming it is E2(r), the channel position E2(r) is determined according to the following formula:

[0013]

[0014] Among them, CH iso is the position of the center of rotation, and R is the number of readings obtained after one rotation of the detector.

[0015] In one embodiment, the step S108 includes: step S117, padding the projection after the convolution filtering with zeros in the frequency domain.

[0016] In one embodiment, the step S108 includes: step S121 , fitting the channel position C(r) corresponding to the center of the phantom.

[0017] In one embodiment, the step S114 includes: step S122, constructing a function ΔF(k, l)=min(∑ r (diff(k,l,r)-a(k,l)·P sim (k, r)-c(k, l)·(P sim (k+1, r)-P sim (k-1, r))) 2 ); Step S124, solve a(k, l) and c(k, l); Step S126, calculate the channel spectrum response difference correction table t according to the following formula CCR (k, l) and channel spacing correction table t spacing (k, l):

[0018] t CCR (k, l) = HP(a(k, l))

[0019] t spacing (k, l) = HP(c(k, l))

[0020] HP() is a high-pass filter function, and filtering is performed in the k direction.

[0021] According to a second aspect of the present invention, a computer program is provided, which can implement the steps of the above method when executed by a processor.

[0022] According to a third aspect of the present invention, there is provided a computer-readable storage medium having a computer program stored thereon, which can implement the steps of the above-mentioned method when executed by a processor.

[0023] According to a fourth aspect of the present invention, a computer tomography apparatus is provided, comprising the computer-readable storage medium as described above.

[0024] The present invention provides a method, computer program, computer-readable storage medium, and computed tomography apparatus for generating a channel spectral response difference correction table and a channel spacing correction table, which are unaffected by residual beam hardening and effectively generate ideal projections. Compared to the prior art, the present invention provides a more accurate table. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] The preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings, so that those skilled in the art will understand the above and other features and advantages of the present invention more clearly. In the accompanying drawings:

[0026] Figure 1The figure is a schematic flow chart of a method for generating a channel spectrum response difference correction table and a channel spacing correction table according to an embodiment of the present invention.

[0027] Figure 2A for Figure 1 Schematic diagram of the original projection of the method shown.

[0028] Figure 2B for Figure 1 Schematic diagram of the average projection method shown.

[0029] Figure 2C for Figure 1 Schematic diagram of the parallel projection method shown.

[0030] Figure 2D for Figure 2C The schematic diagram of the parallel projection after edge enhancement is shown.

[0031] Figure 2E for Figure 1 Schematic diagram of the simulated parallel projection of the method shown.

[0032] Figure 2F for Figure 1 Schematic diagram of the simulated raw projection of the method shown.

[0033] Figure 2G for Figure 1 Schematic diagram of the difference between the original projection and the simulated ideal projection of the method shown.

[0034] Figure 3 Schematic diagram of original projection and filtered projection according to an embodiment of the present invention.

[0035] Figure 4 Schematic diagram of a filtered projection of a partially missing original projection according to an embodiment of the present invention.

[0036] Figure 5 Schematic diagram of the distance between the position of the phantom and the rotation center and its fitting according to an embodiment of the present invention.

[0037] Figure 6 Schematic diagram of the X-ray path in a cylindrical water phantom according to an embodiment of the present invention.

[0038] Figure 7A Schematic diagram of projection difference according to an embodiment of the present invention.

[0039] Figure 7B Schematic diagram of the projection difference calculated by the fitting method in the prior art.

[0040] Figure 7C yes Figure 7A and Figure 7BSchematic diagram of the projection difference signal at a specific reading.

[0041] Figure 8A It is a schematic diagram of a channel spectrum response difference correction table calculated by the fitting method in the prior art.

[0042] Figure 8B FIG. 4 is a schematic diagram of a channel spectrum response difference correction table according to a method for generating a channel spectrum response difference correction table and a channel spacing correction table according to an embodiment of the present invention.

[0043] Figure 9A Schematic diagram of a channel spacing correction table calculated using a fitting method in the prior art.

[0044] Figure 9B 4 is a schematic diagram of a channel spacing correction table according to a method for generating a channel spectrum response difference correction table and a channel spacing correction table according to an embodiment of the present invention. DETAILED DESCRIPTION

[0045] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention is further described in detail with reference to the following examples.

[0046] Figure 1 FIG. 1 is a schematic flow chart of a method 100 for generating a channel spectrum response difference correction table and a channel spacing correction table according to an embodiment of the present invention. Figure 1 As shown, method 100 includes step S102 , step S104 , step S106 , step S108 , step S110 , step S112 , and step S114 .

[0047] In step S102, a raw projection P(k, l, r) of a water phantom is obtained. The water phantom covers all channels of the detector, where k represents the channel number, l represents the sequence number, and r represents the reading number. In this embodiment, a 20 cm water phantom is used, eccentrically positioned relative to the ISO center (e.g., 150 mm apart for a 500 mm FOV). The phantom is positioned radially with a certain tolerance, e.g., 5 mm. The raw projection P(k, l, r) is corrected for air, beam hardening (WBC), and scatter. Figure 2A for Figure 1 A schematic diagram of the original projection P(k, l, r) of method 100 is shown.

[0048] Step S104: average the original projections P(k, l, r) in the row direction to obtain the average projection P mean (k, r). Figure 2B for Figure 1 The average projection P of the method 100 shown mean Schematic diagram of (k, r).

[0049] Step S106: for each channel, interpolate the average projection P in the projection angle direction. mean (k, r), to obtain the parallel projection P paral (k, r). Parallel projection P paral (k, r) are non-isometric. Figure 2C for Figure 1 The parallel projection P of the method 100 shown paral Schematic diagram of (k, r).

[0050] Step S108, according to the parallel projection P paral (k, r), determine the channel position corresponding to the center of the phantom. Figure 1 As shown, in this embodiment, step S108 includes step S116, step S118, and step S120.

[0051] Step S116, parallel projection P paral (k, r) is convolution filtered. For example, convolution filtering can be performed according to the following formula:

[0052] FP=iFFT(FFT(p paral )×H(ω))

[0053] Where FFT() is the Fourier transform, H(ω) is the convolution filter function, and iFFT() is the inverse Fourier transform. H(ω) is the product of the ramp function and the windowing function.

[0054] In order to increase the spatial sampling density, step S108 may further include:

[0055] Step S117: padding the projection after convolution filtering with zeros in the frequency domain.

[0056] After filtering, the edges of the phantom are enhanced, while other parts are weakened. Figure 2D for Figure 2C The parallel projection P shown paral (k, r) Schematic diagram after edge enhancement. Figure 3 Schematic diagram of an original projection P and a filtered projection FP according to an embodiment of the present invention.

[0057] Step S118, determining the channel positions E1(r) and E2(r) corresponding to the edges of the phantom in the filtered parallel projection. Figure 3 As shown, the channel positions E1(r) and E2(r) corresponding to the edge of the phantom are at the local minimum point of the FP curve. paral (k, r) is missing some data, such as Figure 2D As shown, we can only find the channel corresponding to one edge directly. Figure 4Schematic diagram of a filtered projection FP of a partially missing original projection according to an embodiment of the present invention. Assuming that the channel position E1(r) is found, the channel position E2(r) can be determined according to the following formula:

[0058]

[0059] Among them, CH iso is the position of the center of rotation, and R is the number of readings obtained after one rotation of the detector.

[0060] Step S120 , calculating the channel position C(r) corresponding to the center of the phantom based on the channel positions E1(r) and E2(r) corresponding to the edges of the phantom, where C(r)=(E1(r)+E2(r)) / 2.

[0061] In order to improve the accuracy of the phantom position (and the channel position corresponding to the phantom center), step s108 may further include:

[0062] Step S121 , fitting the channel position C(r) corresponding to the center of the phantom. Figure 5 Schematic diagram of the distance between the position of the phantom and the rotation center and its fitting according to an embodiment of the present invention.

[0063] Step S110, calculate the simulated parallel projection according to the channel position corresponding to the center of the phantom The Chinese patent application (application number 202310555264.2) discloses a method for correcting beam hardening projection, a method for determining the position of a phantom, related devices and equipment, including calculating a simulated (theoretical) parallel projection based on the channel position corresponding to the center of the phantom. Figure 6 Schematic diagram of the X-ray path in a cylindrical water phantom according to an embodiment of the present invention. In this embodiment, the phantom is a PMMA water phantom, and the attenuation per unit length of the water phantom can be calculated by the following formula:

[0064] μl=μ1(l1-l2)+μ2l2

[0065] in:

[0066] μ1 and μ2 are the attenuation coefficients of PMMA and water, respectively;

[0067] l1 and l2 are the path lengths of the currently projected X-ray through the water phantom and through the water, respectively;

[0068] μ is the attenuation per unit length of the water model;

[0069] l is the path length through the water model, same as l1.

[0070] According to μ and l, the simulated parallel projection can be calculated

[0071] Step s112: for each channel, interpolate and simulate parallel projection in the projection angle direction To obtain the simulated ideal projection P sim (k, r), and calculate the projection difference diff(k, l, r) between it and the original projection, diff(k, l, r) = P(k, l, r) - P sim (k, r). Figure 7A is a schematic diagram of the projection difference diff(k, l, r) according to an embodiment of the present invention, Figure 7B This is a schematic diagram of the projection difference calculated by the fitting method of the prior art. Figure 7C yes Figure 7A and Figure 7B The projection difference of the signal at specific readings 1 and 2 is shown in Figure 2. The fitting result is affected by the residual WBC, which causes the position of the model to deviate from the actual position. This deviation will be exacerbated when some data are missing. Figure 7A and Figure 7B As shown, readings 1 and 2 are two positions symmetrical about ISO, and the distance between the phantom and the focus is also the same. This means that at these two positions, the measured data and the simulated data should have a profile that is almost symmetrical about ISO. Due to the design of the CT detector's partial fan beam, some data is missing. Figure 7C As shown, at two positions, the fitting results do not have a symmetrical shape, while the results of the method of the present invention do, which means that the fitting method has obvious fitting deviations.

[0072] Step S114: According to the projection difference diff(k, l, r) and the simulated ideal projection P sim (k, r), calculate the channel spectrum response difference correction table t CCR (k, l) and channel spacing correction table t spacing (k, l). Step S114 may include:

[0073] Step S122, construct the function ΔF(k, l)=min(∑ r (diff(k,l,r)-a(k,l)·P sim (k, r)-c(k, l)·(P sim (k+1, r)-P sim (k-1, r))) 2 );

[0074] In step S124, a(k, l) and c(k, l) are solved. a(k, l) and c(k, l) are the channel spectrum response difference correction coefficients and channel spacing correction coefficients, respectively. They can be calculated by minimizing ΔF(k, l). That is:

[0075]

[0076]

[0077] The final linear equation can be obtained from the above equation:

[0078]

[0079]

[0080] By solving the above linear equation, we can calculate a(k, l) and c(k, l).

[0081] Step S126: Calculate the channel spectrum response difference correction table t according to the following formula: CCR (k, l) and channel spacing correction table t spacing (k, l):

[0082] t CCR (k, l) = HP(a(k, l))

[0083] t spacing (k, l) = HP(c(k, l))

[0084] Where HP() is a high-pass function. Since low-frequency offsets do not affect image quality, a(k, l) and c(k, l) are filtered in the k direction using a high-pass filter to obtain the final CCR table and spacing table. Figure 8A is a schematic diagram of a channel spectrum response difference correction table calculated by the fitting method of the prior art, Figure 8B FIG. 4 is a schematic diagram of a channel spectrum response difference correction table according to a method for generating a channel spectrum response difference correction table and a channel spacing correction table according to an embodiment of the present invention. Figure 9A is a schematic diagram of a channel spacing correction table calculated using a fitting method in the prior art. Figure 9B 4 is a schematic diagram of a channel spacing correction table according to a method for generating a channel spectrum response difference correction table and a channel spacing correction table according to an embodiment of the present invention.

[0085] According to a second aspect of the present invention, a computer program is provided. When executed by a processor, the program can implement the steps of the method 100 for generating a channel spectrum response difference correction table and a channel spacing correction table.

[0086] According to a third aspect of the present invention, a computer-readable storage medium is provided, on which a computer program is stored. When the program is executed by a processor, the steps of the method 100 for generating a channel spectrum response difference correction table and a channel spacing correction table can be implemented.

[0087] According to a fourth aspect of the present invention, a computer tomography apparatus is provided, comprising the computer-readable storage medium as described above.

[0088] The present invention provides a method, computer program, computer-readable storage medium, and computed tomography apparatus for generating a channel spectral response difference correction table and a channel spacing correction table, which are unaffected by residual beam hardening and effectively generate ideal projections. Compared to the prior art, the present invention provides a more accurate table.

[0089] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for generating a channel spectrum response difference correction table and a channel spacing correction table, comprising: Step S102: obtaining an original projection P(k, l, r) of a water phantom, where the water phantom covers all channels of the detector, where k represents the channel number, l represents the sequence number, and r represents the reading number; Step S104: average the original projections P(k, l, r) in the row direction to obtain the average projection P mean (k,r); Step S106: for each channel, interpolate the average projection P in the projection angle direction. mean (k,r), to obtain the parallel projection P paral (k,r); Step S108, according to the parallel projection P paral (k, r), determine the channel position corresponding to the center of the phantom; Step S110, calculating the simulated parallel projection according to the channel position corresponding to the center of the phantom Step S112: for each channel, interpolate and simulate parallel projection in the projection angle direction To obtain the simulated ideal projection P sim (k,r), and calculate the projection difference diff(k,l,r) between it and the original projection, diff(k,l,r)=P(k,l,r)-P sim (k,r); Step S114: According to the projection difference diff(k,l,r) and the simulated ideal projection P sim (k, r), calculate the channel spectrum response difference correction table t CCR (k, l) and channel spacing correction table t spacing (k,l).

2. The method according to claim 1, wherein: The step S108 includes: Step S116, parallel projection P paral (k,r) performs convolution filtering; Step S118, determining channel positions E1(r) and E2(r) corresponding to the edges of the phantom in the filtered parallel projections; Step S120 , calculating the channel position C(r) corresponding to the center of the phantom based on the channel positions E1(r) and E2(r) corresponding to the edges of the phantom, where C(r)=(E1(r)+E2(r)) / 2.

3. The method according to claim 2, wherein: The step S116 includes performing convolution filtering according to the following formula: FP=iFFT(FFT(P paral )×H(ω)) Among them, FFT() is Fourier transform, H(ω) is the convolution filter function, and iFFT() is the inverse Fourier transform.

4. The method according to claim 3, wherein: H(ω) is the product of the ramp function and the windowing function.

5. The method according to claim 2, wherein: The parallel projection P paral The missing data of (k, r) results in the loss of one of the channel positions E1(r) and E2(r) corresponding to the edge of the phantom. Assuming it is E2(r), the channel position E2(r) is determined according to the following formula: Among them, CH iso is the position of the center of rotation, and R is the number of readings obtained after one rotation of the detector.

6. The method according to claim 2, wherein: The step S108 includes: Step S117: padding the projection after convolution filtering with zeros in the frequency domain.

7. The method according to claim 2, wherein: The step S108 includes: Step S121 , fitting the channel position C(r) corresponding to the center of the phantom.

8. The method according to claim 1, wherein: The step S114 includes: Step S122, construct the function ΔF(k,l)=min(∑ r (diff(k,l,r)-a(k,l)·P sim (k,r)-c(k,l)·(P sim (k+1,r)-p sim (k-1,r))) 2 ); Step S124, solving a(k, l) and c(k, l); Step S126: Calculate the channel spectrum response difference correction table t according to the following formula: CCR (k, l) and channel spacing correction table t spacing (k,l): t CCR (k,l)=HP(a(k,l)) t spacing (k,l)=HP(c(k,l)) HP() is a high-pass filter function, and filtering is performed in the k direction.

9. A computer program, which, when executed by a processor, can implement the steps of the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the steps of the method according to any one of claims 1 to 8 can be implemented.

11. A computer tomography apparatus comprising the computer-readable storage medium according to claim 10.

Images