Method for generating channel spectral response difference correction table and channel pitch correction table, computer program

Abstract

The invention discloses a method for generating a channel spectral response difference correction table and a channel spacing correction table, a computer program, a computer readable storage medium and computed tomography equipment. The method for generating the channel spectral response difference correction table and the channel spacing correction table comprises the following steps: determining a channel position corresponding to the center of a motif by using parallel projection Pparal (k, r), calculating simulation parallel projection # imgabs0 # simulation ideal projection Psim (k, r) according to the channel position, and finally calculating the channel spectral response difference correction table tCCR (k, l) and the channel spacing correction table tspacing (k, l) according to the projection difference diff (k, l, r) and the simulation ideal projection Psim (k, r). According to the method for generating the channel spectral response difference correction table and the channel spacing correction table, the computer program, the computer readable storage medium and the computed tomography equipment, the influence of beam hardening residues is avoided, and the ideal projection can be effectively generated.

Description

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

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

[0003] The X-ray beam generated by the X-ray tube covers a wide spectrum that changes as the attenuation of the object increases. The different spectral responses of each detector channel will cause signal offset, resulting in ring artifacts on the image. In addition, due to mechanical errors, the inevitable non-equidistant detector units may form ring and strip artifacts on the image. To prevent these artifacts, channel spectral response difference correction (CCR) and channel spacing correction (or Spacing) are required.

[0004] The key to generating the channel spectrum 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, but due to the inevitable beam hardening residue, the fitting method makes the fitting result deviate from the actual ideal projection, especially when part of the projection is missing. Summary of the invention

[0005] In view of this, the present invention proposes 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.

[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 model, wherein the water model can cover all channels of the detector, wherein k represents the channel number, l represents the sorting number, and r represents the reading number; step S104, averaging the original projection 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 a parallel projection operation on 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 model in the parallel projection after filtering; step S120, calculating the channel position C(r) corresponding to the center of the model according to the channel positions E1(r) and E2(r) corresponding to the edges of the model, 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 missing of one of the channel positions E1(r) and E2(r) corresponding to the edge of the model. 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, solving a(k, l) and c(k, l); Step S126, calculating 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, there is provided a computer tomography device comprising the computer readable storage medium as described above.

[0024] 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 beam hardening residues and can effectively generate ideal projections. Compared with the prior art, the method of the present invention can provide 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 the above and other features and advantages of the present invention will be more apparent to those skilled in the art. In the accompanying drawings:

[0026] Figure 1The present invention is a schematic flowchart 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 representation of the average projection of the 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] Fig. 7A Schematic diagram of projection difference according to an embodiment of the present invention.

[0039] Figure 7B This is a schematic diagram of the projection difference calculated by the fitting method in the prior art.

[0040] Figure 7C yes Fig. 7A and Figure 7BSchematic representation of the projection of the difference signal at a particular reading.

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

[0042] Figure 8B 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] Fig.9A It is a schematic diagram of a channel spacing correction table calculated by a fitting method in the prior art.

[0044] Fig. 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 embodiments.

[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] Step S102, obtain an original projection P(k, l, r) of a water phantom, which can cover all channels of the detector, where k represents the channel number, l represents the sorting number, and r represents the reading number. This embodiment uses a 20 cm water phantom, which is placed eccentrically relative to the rotation center (ISOcenter), for example, 150 mm apart under a FOV of 500 mm. The water phantom is placed radially with a certain tolerance, for example, 5 mm. The original projection P(k, l, r) is corrected for air, beam hardening correction (WBC) based on the water phantom, and scattering. Figure 2A for Figure 1 A schematic diagram of the original projection P(k, l, r) of method 100 is shown.

[0048] Step S104, averaging the original projections P(k, l, r) in the row direction to obtain an 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 model. 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 convolutionally filtered. For example, convolution filtering can be performed according to the following formula:

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

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

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

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

[0056] After filtering, the edges of the model are enhanced and 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 edges of the model are at the local minimum points 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 according to 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 center of the phantom), step s108 may also 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, calculating 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, and related devices and equipment, including calculating and simulating (theoretical) parallel projections 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 model 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, we can calculate the simulated parallel projection

[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). Fig. 7A is a schematic diagram of a projection difference diff(k, l, r) according to an embodiment of the present invention, Figure 7B is a schematic diagram of the projection difference calculated by the fitting method of the prior art, Figure 7C yes Fig. 7A and Figure 7B Schematic diagram of the projection difference of the signal at specific reads 1 and 2. The fitting result is affected by the WBC carryover, causing the motif position to deviate from the actual position, and this deviation will worsen when some data are missing. Fig. 7A and Figure 7B As shown in the figure, 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 an obvious fitting deviation.

[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, 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 );

[0074] Step S124, solve a(k, l) and c(k, l), a(k, l) and c(k, l) are channel spectrum response difference correction and channel spacing correction coefficients respectively, which 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] 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. Fig. 8A is a schematic diagram of a channel spectrum response difference correction table calculated by a fitting method of the prior art, Figure 8B 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. Fig.9A is a schematic diagram of a channel spacing correction table calculated by a fitting method of the prior art, Fig. 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, which, when executed by a processor, 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, which, when executed by a processor, can implement the steps of the method 100 for generating a channel spectrum response difference correction table and a channel spacing correction table.

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

[0088] 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 beam hardening residues and can effectively generate ideal projections. Compared with the prior art, the method of the present invention can provide 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 principle of the present invention should be included in the protection scope 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 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: average 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, 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, characterized in that: The step S108 comprises: 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 according to the channel positions E1(r) and E2(r) corresponding to the edges of the phantom, C(r)=(E1(r)+E2(r)) / 2.

3. The method according to claim 2, characterized in that: 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, characterized in that: H(ω) is the product of the ramp function and the windowing function.

5. The method according to claim 2, characterized in that: The parallel projection P paral The missing data of (k, r) results in the missing of one of the channel positions E1(r) and E2(r) corresponding to the edge of the model. 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, characterized in that: The step S108 includes: Step S117, in the frequency domain, padding the projection after the convolution filtering with zeros.

7. The method according to claim 2, characterized in that: 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, characterized in that: The step S114 comprises: 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 described in 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