CBCT large-view cone beam artifact elimination method and device

Abstract

The invention relates to a CBCT large-view cone beam artifact elimination method. The method comprises the steps that S1, cone beam CT projection data are collected and corrected; s2, performing parallel beam rearrangement on the projection data according to the structure of the cone beam scanning system; s3, determining a z-direction expansion reconstruction range according to the reconstruction view and the geometric structure; s4, designing a weighting function, and performing smooth transition between the full scanning area and the partial scanning area; and S5, through single convolution and back projection, generating a taper angle artifact removal image for expanding the z range. According to the method, the z-axis scanning range is remarkably expanded by reconstructing the 180-degree coverage voxels at the top and the bottom in the Z direction, projection data are fully utilized, and the ineffective radiation dosage is reduced; real-time reconstruction is realized through single convolution design and GPU parallel acceleration, reconstruction calculation complexity is reduced, image discontinuity is reduced through smooth transition weight, and diagnosis reliability is improved; the core bottleneck of traditional cone beam CT in expansion reconstruction is solved, and a high-coverage, low-dose and high-quality imaging solution is provided for clinic.

Description

Technical Field

[0001] The present invention relates to the field of oral imaging technology, and particularly to a method and device for eliminating large-field cone-beam artifacts in CBCT. Background Art

[0002] Oral CBCT (Cone Beam Computed Tomography) is a medical technology specifically used for oral and maxillofacial imaging, which can provide high-resolution three-dimensional images to help dentists more accurately diagnose diseases of teeth, jaws and related structures. The circular trajectory scanning method has become one of the most widely used scanning modes in CBCT due to its simple engineering implementation and convenient control. With the progress of flat panel detector technology, the detector size has been continuously increasing, and the scanning range has been significantly expanded. However, the filtered back-projection reconstruction algorithm is accurate for reconstruction when the cone angle is zero, and when the cone angle increases, the algorithm switches to approximate reconstruction, resulting in an increase in reconstruction error and serious cone-beam artifacts appearing at the top and bottom of the image. This problem is particularly prominent in large cone angle scanning, mainly manifested as longitudinal intensity attenuation, which affects the diagnostic accuracy.

[0003] Traditional two-dimensional X-ray imaging has limitations in oral diagnosis, such as low resolution and lack of three-dimensional information. CBCT, through the combination of a conical X-ray beam and a flat panel detector, can quickly obtain high-quality three-dimensional images of complex structures at low radiation doses. Although the FDK algorithm has become the mainstream reconstruction method due to its high efficiency, its approximate characteristics lead to image deviation, especially when the detector size increases and the cone angle increases, the artifact problem becomes more significant. How to effectively eliminate the artifacts generated by large cone angle reconstruction has become a key challenge in improving the imaging quality of oral CBCT.

[0004] For the cone-beam artifacts generated in the circular orbit cone-beam flat panel CT imaging of large-size detector arrays, the existing solutions can be divided into the following four categories:

[0005] 1. 3D Weighting Function Method

[0006] This method eliminates the cone-beam artifacts caused by large cone angles by designing different 3D weighting functions. Specifically, it can be divided into two implementation methods: pre-filtering weighting, where the weighting function is based on the fan angle position and detector distance of the reconstruction point, but due to the dynamic change of the weights in some scanning regions among voxels, the computational complexity increases significantly. Post-filtering weighting, where the weighting function is applied after data filtering, reducing the computational complexity and separately processing full scans and local scans, but it is not applicable to short-scan CT imaging.

[0007] 2. Data Rearrangement Method

[0008] The data rearrangement method, based on the circular orbit cone-beam FDK algorithm, rearranges the fan beam into a parallel beam through bilinear interpolation, changes the filtering path to compensate for the longitudinal intensity loss of the image, thereby improving the image quality.

[0009] Although the data rearrangement method can compensate for the loss at the top of the cone angle, it will change the filtering direction, increase the calculation time, and may introduce noise, thus affecting the image quality.

[0010] 3. Iterative reconstruction algorithm

[0011] The iterative reconstruction algorithm is particularly suitable for scenarios with insufficient projection data, missing projection angles, or uneven projection intervals.

[0012] However, the iterative reconstruction algorithm has a slow reconstruction speed. Its time consumption is closely related to the number of projections, the size of the reconstructed image, and the number of iterations, making it difficult to meet the clinical real-time requirements.

[0013] 4. Deep learning method

[0014] The deep learning-based method can automatically learn features from the training data. Without the need for manual design of complex feature extraction algorithms, it can effectively eliminate cone-beam artifacts.

[0015] However, the deep learning method requires a large amount of accurate and reliable training data, and the training process is complex, time-consuming, and requires high computing resources, which may affect the model performance. Summary of the Invention

[0016] Aiming at the deficiencies of the existing methods for eliminating large-field-of-view cone-beam artifacts in CBCT, the present invention designs a CBCT large-field-of-view cone-beam artifact elimination device and method. The reconstruction range of the traditional FDK algorithm is limited by the size of the flat panel detector and can only reconstruct the area covered by 360°, resulting in limited longitudinal (z-axis) range and low dose utilization. In addition, the existing cone angle artifact elimination and reconstruction methods have problems such as high computational complexity, and unclear elimination of noise and artifacts. The present invention aims to overcome the limitations of the existing technology and provide a CBCT large-field-of-view cone angle artifact elimination method with high computational efficiency, wide applicability, strong real-time performance, and low resource requirements to meet the actual clinical needs.

[0017] To solve the above technical problems, the present invention provides a CBCT large-field-of-view cone-beam artifact elimination method, including the following steps:

[0018] Step S1: Collect cone-beam CT projection data and perform correction;

[0019] Step S2: Rearrange the projection data into parallel beams according to the cone-beam scanning system structure;

[0020] Step S3: Determine the z-direction extended reconstruction range according to the reconstruction field of view and geometric structure;

[0021] Step S4: Design a weight function to achieve smooth transition between the full scan area and the partial scan area;

[0022] Step S5: Generate an image with cone angle artifacts removed and an extended z range through single convolution and backprojection.

[0023] Further, in step S1, the original projection data is collected by rotating the CBCT device 360°, and physical correction is performed on the original projection data.

[0024] Further, in step S1, the correction of the original projection data includes dark field correction, bad pixel correction, detector response correction, and scatter correction, and geometric correction is performed according to the system geometric parameters.

[0025] Further, in step S2, in the circular trajectory cone beam scanning system, the light source rotates around the rotation axis (i.e., the z-axis) with a radius R sid to rearrange the cone beam projection data into parallel data.

[0026] Further, in step S2, during the data rearrangement process, several adjacent points are selected from the cone beam projection data for interpolation using the Lagrange interpolation method to approximately obtain the parallel projection data.

[0027] Further, in step S3, by determining the distance R from the ray source to the rotation center sid and the radius R of the reconstructed axial plane, the reconstruction range in the z direction is extended by times.

[0028] Further, in step S4, the formula for the weight function ω(r, β) is:

[0029] ω(r, β) = (1 - ω r (r, β)) * ω f (r, β) + ω r (r, β) * ω β (r, β);

[0030] Where:

[0031] ω β (r, β) is the partial scan weight function, determined according to the projection angle β and the reconstructed voxel position r = (x, y, z);

[0032] ω r (r, β) is the transition weight function, used to smoothly transition between the full scan area and the partial scan area;

[0033] ω f (r, β) is the full scan weight function, ω f (r, β) = 1;

[0034] The full-scan weight function is directly used in the central region (full-scan region), gradually switched to the partial-scan weight function in the transition region, and the weights are dynamically adjusted in the extended region (partial-scan region) to compensate for the insufficient data redundancy.

[0035] Furthermore, the partial-scan weight function ω β (r, β) is calculated by the formula:

[0036]

[0037] where is a smooth transition function, β start is the starting projection angle, and Δβ is the projection angle correction amount.

[0038] Furthermore, the transition weight function ω r (r, β) is calculated by the formula:

[0039]

[0040] where R sid is the distance from the ray source to the rotation center, represents the cone angle, z represents the reconstruction field of view size in the z direction, r represents the distance from the reconstructed voxel to the rotation center, is a smooth transition function.

[0041] Furthermore, in step S5, single-filter back-projection reconstruction is performed on the rearranged parallel projection data in combination with the weight function ω(r, β), and the specific formula is as follows:

[0042]

[0043] where P(β, t, s) is the rearranged projection data, k(ξ) is the filter kernel, and f(x, y, z) is the reconstructed voxel.

[0044] The present invention also provides a CBCT large-field-of-view cone-beam artifact elimination device, including:

[0045] at least one processor; and

[0046] at least one memory communicatively connected to the processor;

[0047] wherein the memory stores instructions executable by the processor, and the instructions are executed by the processor to enable the device to perform the aforementioned CBCT large-field-of-view cone-beam artifact elimination method.

[0048] Advantages of the present invention:

[0049] Compared with the prior art, the CBCT large-field-of-view cone-beam artifact elimination method of the present invention has the following technical effects:

[0050] (1) For large flat-panel detectors, the z-direction reconstruction field of view of the traditional FDK reconstruction method is limited. The present invention can significantly increase the longitudinal field of view through z-direction weighted expansion, providing a larger field-of-view imaging for orthodontists in three-dimensional orthodontics.

[0051] (2) For partial scan data at the top and bottom of the detector, the present invention can make full use of 180° scan data, achieve 100% dose utilization, and reduce the CTDI value (CT radiation dose).

[0052] (3) Compared with other existing cone-angle artifact elimination reconstruction methods that have problems such as unobvious noise and artifact elimination, the present invention can better compensate for the inaccurate CT values caused by large cone angles and the dark areas at the top and bottom in the Z direction through smoothing weights, restore the true HU value (Hounsfield Unit), and can better control the influence of noise.

[0053] (4) The present invention uses dynamic weighting and single-filter backprojection, and through GPU parallel acceleration, greatly reduces the reconstruction calculation complexity, can achieve fast imaging, and does not require a large amount of training data. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] The following further elaborates on the specific embodiments of the present invention with reference to the drawings.

[0055] Figure 1 It is a flowchart of the CBCT large-field-of-view cone-beam artifact elimination method of the present invention.

[0056] Figure 2 It is a schematic diagram of the circular trajectory cone-beam scanning geometry.

[0057] Figure 3 It is a coordinate schematic diagram of cone-beam data P1(θ, m, n) and rearranged data P(β, t, s).

[0058] Figure 4 It is a system structure diagram of the axial plane and the sagittal plane.

[0059] Figure 5 It is a schematic diagram of the distribution of the cone-beam reconstruction area. DETAILED DESCRIPTION OF THE INVENTION

[0060] Example 1

[0061] Combined with Figure 1 , the CBCT large-field-of-view cone-beam artifact elimination method of this embodiment specifically includes the following steps:

[0062] Step S1: Collect cone-beam CT projection data and perform correction.

[0063] Specifically, in this embodiment, in step S1, the original projection data is collected by rotating the CBCT device by 360°, and physical correction is performed on the original projection data.

[0064] Preferably, in this embodiment, in step S1, the correction of the original projection data includes dark field correction, bad pixel correction, detector response correction, and scatter correction, and geometric correction is performed according to the system geometric parameters to ensure the accuracy of the original data.

[0065] Step S2: Rearrange the projection data into parallel beam data according to the structure of the cone beam scanning system.

[0066] Preferably, in this embodiment, as Figure 2 shown, in step S2, in the circular trajectory cone beam scanning system, the light source rotates around the rotation axis (i.e., the z-axis) with a radius R sid to rotate, and the cone beam projection data P1(θ, m, n) is rearranged into parallel data P(β, t, s). The rearrangement formula is as follows:

[0067]

[0068] where θ represents the projection angle, (m, n) represents the coordinates on the detector, β represents the projection angle after rearrangement, (t, s) represents the detector coordinates of the projection data after rearrangement, and R sid represents the distance from the ray source center to the rotation center.

[0069] Specifically, in this embodiment, as Figure 3 shown, during the process of rearranging the cone beam projection data P1(θ, m, n) into parallel projection data P(β, t, s), it is impossible to ensure that each point of the parallel projection data can be found in the cone beam projection data. Therefore, during the data rearrangement process, the Lagrange interpolation method needs to be used to approximately obtain the parallel projection data.

[0070] Preferably, in this embodiment, in step S2, during the data rearrangement process, several adjacent points of the cone beam projection data are selected for interpolation by the Lagrange interpolation method to approximately obtain the parallel projection data.

[0071] Specifically, in this embodiment, the Lagrange interpolation method is a method of constructing a polynomial through known points. Given n different points (x0, y0), (x1, y1) …… (x n-1 , y n-1 ), a polynomial P(x) with a degree not exceeding n - 1 can be found such that P(x i ) = y i holds for all i.

[0072] In this embodiment, for the cone-beam projection data P1(θ, m, n), five adjacent points are selected for interpolation. First, the Lagrange basis functions are constructed, and the specific formulas are as follows:

[0073]

[0074] Specifically, each term of the Lagrange basis function is multiplied by the corresponding y i and accumulated, and the specific formula is as follows:

[0075]

[0076] Step S3: Determine the z-direction extended reconstruction range according to the reconstruction field of view and geometric structure. Specifically, in this embodiment, as Figure 4 shown, the longitudinal range of the detector is limited, and there are some slices far from the central plane that are not fully irradiated. When z is large enough, some positions (x, y) in the slice cannot be measured. The irradiated voxel region and the non-irradiated region are divided by a boundary line. The slice region close to the focal point is not irradiated, while the voxels at the far end of the boundary line are irradiated. Based on symmetry, the case of z > 0 is taken as an example in this embodiment. Specifically, the reconstruction boundary conditions for the full-scan region (360° scan) and the partial-scan region (scan greater than 180° and less than 360°) are defined, and the specific formulas are as follows:

[0077]

[0078] where, R sid is the distance from the ray source to the rotation center, β represents the projection angle, represents the cone angle, and z represents the size of the reconstruction field of view in the z direction.

[0079] Specifically, according to Figure 4 it can be seen that the condition for the full-scan region (360° scan) is:

[0080]

[0081] The condition for the partial-scan region (scan greater than 180° and less than 360°) is:

[0082]

[0083] According to the formula transformation, the highest increment of the z-direction extended reconstruction can be determined:

[0084]

[0085] Preferably, in this embodiment, in step S3, by determining the distance R sid from the ray source to the rotation center and the radius R of the reconstructed axial plane, the z-direction extended reconstruction range is increased by times.

[0086] Step S4: Design a weight function to achieve a smooth transition between the full-scan region and the partial-scan region.

[0087] Specifically, in this embodiment, as Figure 5 shown, region a is the full-scan region (360° scan) and can be fully reconstructed; region b is the partial-scan region (scan angle greater than 180° and less than 360°) and can be reconstructed by weighted patching with weights; region c is the region with a scan angle less than 180°, and the weight function cannot perform patching reconstruction.

[0088] Specifically, in this embodiment, the partial-scan weight function ω β (r, β) is determined according to the projection angle β and the reconstructed voxel position r = (x, y, z), and its specific calculation formula is:

[0089]

[0090] where is a smooth transition function, β start is the starting projection angle, and Δβ is the projection angle correction amount.

[0091] Specifically, in this embodiment, the projection angle correction amount Δβ is determined by the voxel position and the scanning geometry, and the specific formula is as follows:

[0092]

[0093] where R sid is the distance from the ray source to the rotation center, represents the cone angle, z represents the reconstruction field of view size in the z direction, and r represents the distance from the reconstructed voxel to the rotation center.

[0094] Specifically, in this embodiment, to avoid sudden changes between the full-scan and partial-scan regions, a transition weight function ω r (r, β) is introduced, and the calculation formula is as follows:

[0095]

[0096] where R sid is the distance from the ray source to the rotation center, represents the cone angle, z represents the reconstruction field of view size in the z direction, r represents the distance from the reconstructed voxel to the rotation center, is a smooth transition function.

[0097] Specifically, in this embodiment, the full-scan weight and the partial-scan weight are smoothly transitioned within the r±Δr region; it is ensured that the weight changes continuously in the extended region (when z is large), avoiding sudden changes in noise and artifacts. Finally, considering the influence of factors such as the projection angle and the position of the reconstructed voxel, through the transition weight function ω r (r, β) mixes the partial-scan weight ω β (r, β) and the full-scan weight ω f (r, β) = 1, and the calculation formula of the final weight function ω(r, β) is:

[0098] ω(r, β) = (1 - ω r (r, β)) * ω f (r, β) + ω r (r, β) * ω β (r, β);

[0099] In the central region (full-scan region), the full-scan weight function is directly used. In the transition region, it is gradually switched to the partial-scan weight function. In the extended region (partial-scan region), the weight is dynamically adjusted to compensate for the insufficient data redundancy.

[0100] Step S5: Generate a cone-beam artifact-removed image with an extended z range through single convolution and backprojection.

[0101] Preferably, in this embodiment, in step S5, single-filtered backprojection reconstruction is performed on the rearranged parallel projection data in combination with the weight function ω(r, β), and the specific formula is as follows:

[0102]

[0103] where P(β, t, s) is the rearranged projection data, k(ξ) is the filter kernel, and f(x, y, z) is the reconstructed voxel.

[0104] The CBCT large-field-of-view cone-beam artifact elimination method of this embodiment breaks through the longitudinal field-of-view limitation of traditional reconstruction by developing an adaptive weighted extension algorithm in the z-axis direction. Through experimental verification, this technology can expand the effective imaging field-of-view by about 40%, providing the ability to perform a complete three-dimensional reconstruction of the dental arch in the field of orthodontics and supporting clinicians in obtaining more comprehensive occlusal relationship data.

[0105] The CBCT large-field-of-view cone-beam artifact elimination method of this embodiment innovatively realizes the full-dose utilization technology of 180° scan data. Through the detector edge data compensation algorithm, the invalid data regions at the top / bottom are reduced to less than 5%.

[0106] The CBCT large-field-of-view cone-beam artifact elimination method of this embodiment constructs a cone angle correction model based on multi-scale smoothing weights, effectively solving the problem of CT value offset caused by large cone angle scanning. It eliminates the dark area artifacts at the end of the z-axis, and improves the gray scale uniformity to more than 95%.

[0107] The CBCT large-field-of-view cone-beam artifact elimination method of this embodiment can quickly achieve cone angle artifact elimination through dynamic weight allocation and single-filter backprojection, combined with a GPU heterogeneous computing framework, and does not require the support of training data for deep learning models.

[0108] In summary, the CBCT large-field-of-view cone-beam artifact elimination method of the present invention significantly expands the z-axis scanning range by reconstructing 180° coverage voxels at the top and bottom in the Z direction, makes full use of projection data, and reduces the ineffective radiation dose. The single-convolution design and real-time reconstruction are achieved through GPU parallel acceleration, greatly reducing the reconstruction calculation complexity. The smoothly-transitioning weights reduce image discontinuity, greatly improving the diagnostic reliability. Through the innovative weight strategy and efficient calculation process, the present invention solves the core bottleneck in the extended reconstruction of traditional cone-beam CT, providing a high-coverage, low-dose, and high-quality imaging solution for clinical applications.

[0109] Embodiment 2

[0110] This embodiment provides a CBCT large-field-of-view cone-beam artifact elimination device, specifically including:

[0111] At least one processor; and

[0112] At least one memory communicatively connected to the processor;

[0113] Wherein, the memory stores instructions executable by the processor, and when the instructions are executed by the processor, the device is enabled to execute the CBCT large-field-of-view cone-beam artifact elimination method in Embodiment 1.

[0114] Many specific details are set forth in the above description to facilitate a full understanding of the present invention. However, the above description is only a preferred embodiment of the present invention, and the present invention can be implemented in many other ways different from those described herein. Therefore, the present invention is not limited by the specific implementations disclosed above. At the same time, any person skilled in the art can make many possible changes and modifications to the technical solution of the present invention, or modify it into an equivalent embodiment with equivalent changes, without departing from the scope of the technical solution of the present invention. Any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall still fall within the scope of the protection of the technical solution of the present invention.

Claims

1. A method for eliminating cone-beam artifacts in a large field of view of CBCT, characterized in that: It includes the following steps: Step S1: Acquire cone-beam CT projection data and perform correction; Step S2: Rearrange the projection data into parallel beam according to the structure of the cone-beam scanning system; Step S3: Determine the extended reconstruction range in the z direction according to the reconstruction field of view and geometric structure; Step S4: Design a weight function to achieve smooth transition between the full-scan area and the partial-scan area; Step S5: Generate a cone-angle artifact-removed image with an extended z range through single convolution and back-projection.

2. The CBCT large-field cone beam artifact elimination method according to claim 1, wherein: In step S1, the correction of the original projection data includes dark-field correction, bad-pixel correction, detector response correction, and scatter correction, and geometric correction is performed according to the system geometric parameters.

3. The CBCT large field of view cone beam artifact elimination method according to claim 1, wherein: In step S2, in the circular trajectory cone beam scanning system, the light source rotates around the rotation axis with a radius R sid to rearrange the cone beam projection data into parallel data.

4. The CBCT large field of view cone beam artifact elimination method according to claim 3, characterized in that: In step S2, during the data rearrangement process, several adjacent points are selected from the cone-beam projection data and Lagrange interpolation method is used for interpolation to approximately obtain parallel projection data.

5. The CBCT large-field cone beam artifact elimination method according to claim 1, characterized in that: In step S3, by determining the distance R from the ray source to the rotation center sid and the radius R of the reconstructed axial plane, the extended reconstruction range in the z direction is increased times.

6. The CBCT large field of view cone beam artifact elimination method according to claim 1, characterized in that: In step S4, the formula of the weight function ω(r, β) is: ω(r, β) = (1 - ω r (r, β)) * ω f (r, β) + ω r (r, β) * ω β (r, β); Where: ω β (r, β) is a partial scan weight function, which is determined according to the projection angle β and the reconstructed voxel position r; ω r (r, β) is a transition weight function used to smoothly transition between the full scan area and the partial scan area; ω f (r, β) is the full scan weight function; The full-scan weight function is directly used in the central area, gradually switched to the partial-scan weight function in the transition area, and the weight is dynamically adjusted in the extended area to compensate for the insufficient data redundancy.

7. The CBCT large field of view cone beam artifact elimination method according to claim 6, wherein: Partial scan weight function ω β (r, β) is calculated by the formula: Among them, is a smooth transition function, β start is the starting projection angle, and Δβ is the projection angle correction amount.

8. The CBCT large-field-of-view cone beam artifact elimination method according to claim 6, wherein: Transition weight function ω r (r, β) is calculated by the formula: Among them, R sid is the distance from the ray source to the rotation center, represents the cone angle, z represents the reconstruction field of view size in the z direction, r represents the distance from the reconstructed voxel to the rotation center, is a smooth transition function.

9. The CBCT large-field cone-beam artifact elimination method according to claim 1, wherein: In step S5, single-filtered back-projection reconstruction is performed on the rearranged parallel projection data in combination with the weight function ω(r, β), and the specific formula is as follows: f(x, y, z) = ∫0 2π ω(r, β) * P(β, t, s) * k(ξ) dβ; Where P(β, t, s) is the rearranged projection data, k(ξ) is the filter kernel, and f(x, y, z) is the reconstructed voxel.

10. A CBCT large-field cone-beam artifact elimination device, characterized in that: It includes: At least one processor; And At least one memory communicatively connected to the processor; Wherein, the memory stores instructions executable by the processor, and when the instructions are executed by the processor, the device is enabled to execute the CBCT large-field-of-view cone-beam artifact elimination method according to any one of claims 1-9.

Images