CBCT rigid motion artifact correction method and system

Abstract

The invention discloses a CBCT rigid motion artifact correction method and system, and the CBCT rigid motion artifact correction method based on a 3D-2D micro-registration framework outputs a final CT image x and the rigid disturbance quantity of all CT scanning angles through nine steps. According to the method, the estimated reconstructed image is registered to the measurement projection to estimate the rigid disturbance quantity generated in the CBCT scanning process, the rigid disturbance quantity comprises mechanical motion and human body autonomous / non-autonomous motion, the estimated rigid disturbance quantity is used for adjusting imaging geometry so as to reconstruct more accurate images, and the images are also used for achieving more accurate rigid disturbance estimation. The method has the beneficial effects that 1, a rigid motion modeling method is adopted, and a transformation matrix is strictly constrained to meet rigid transformation definition; 2, constructing a gradient-based 3D-2D registration frame, and estimating rigid disturbance by registering a reconstructed image to a measurement projection; and 3, ensuring that the image content is fixed in a specific imaging area by constraining rigid disturbance.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of 3D-2D rigid registration and computed tomography, in particular to a CBCT rigid motion artifact correction method and system based on a 3D-2D differentiable registration framework. BACKGROUND

[0002] CBCT captures volumetric data in a single rotation using a cone-beam of X-rays, providing high-resolution three-dimensional images while minimizing radiation exposure. This characteristic makes CBCT play an important role in many medical fields such as dentistry, maxillofacial surgery, otolaryngology, radiology, and image-guided interventions, providing significant advantages for diagnosis and treatment. Despite its many advantages, CBCT imaging assumes that the subject remains stationary during the scanning process. However, low frame rate flat panel detectors and gantry design limit CBCT scanning speed, increasing the likelihood of subject motion during long scanning processes. Subject rigid motion changes the spatial relationship between the X-ray source, flat panel detector, and the subject, resulting in inaccurate imaging geometry and introducing rigid motion artifacts in the reconstructed images. These artifacts blur anatomical details, making accurate diagnosis difficult and potentially leading to misjudgment of the subject's condition by the doctor. For example, in maxillofacial surgery, rigid motion artifacts can hinder accurate assessment of the severity of fractures or the margins of tumors, affecting surgical planning. In image-guided interventions, rigid motion artifacts compromise the accuracy of surgical navigation, increasing the risk of surgery and reducing the success rate. In addition, the presence of rigid artifacts often requires repeated scanning to obtain images of diagnostic quality, exposing the subject to additional radiation doses. This problem is particularly concerning considering the cumulative effects of radiation exposure and its potential long-term health risks.

[0003] Hardware-based methods aim to reduce subject motion or monitor subject motion through immobilizers. Immobilizers such as headrests, bite blocks, and chin supports limit subject motion during scanning, but cannot completely eliminate it and are not suitable for all subjects, such as stroke patients. Another approach is to use motion tracking systems to provide real-time data of patient motion. However, these systems require additional hardware such as optical or electromagnetic sensors, which increases cost and can complicate the scanning process. In addition, increasing the gantry rotation speed or using multiple X-ray sources can reduce scanning time, effectively "freezing" patient motion and reducing the likelihood of motion. However, this approach is limited by existing hardware and technical capabilities.

[0004] To address the rigid motion artifact problem in CBCT imaging, various software solutions have been proposed, which can be divided into two categories. The first category focuses on data recovery, aiming to correct the projections or images. For example, mirror projection alignment is used to repair the misalignment of projections caused by rigid motion; Fourier consistency of forced chord domain is used to suppress rigid motion artifacts; motion-damaged CT images are modeled as the convolution of a blur kernel and a clean image, and a blind deblurring algorithm is used for image recovery. Inspired by the powerful processing capability of deep learning in medical imaging, researchers also use deep learning to model the mapping from artifact images to clean images, or to learn non-reference indicators to assist image recovery. However, these methods are only post-processing techniques and do not address the root cause of CBCT rigid motion artifacts, so they are difficult to completely correct the anatomical structure deformation caused by them. The second category of methods estimates the motion trajectory of the object during CBCT scanning, and then repairs the imaging geometry for motion compensation reconstruction. These methods mainly include projection domain, image domain and 3D-2D rigid registration. Projection domain optimization methods refine the imaging geometry by maximizing the consistency of the projection data, often using epipolar geometry consistency conditions, Fourier consistency conditions and Helgason-Ludwig consistency conditions. Image domain optimization methods, called autofocus techniques, identify the correct imaging geometry by evaluating the quality of the reconstructed image, often using image quality indicators or deep learning-based non-reference indicators. In addition, 3D-2D rigid registration methods use prior knowledge of rigid motion artifacts to register the reconstructed image to the measured projection data for motion estimation and motion compensation. However, most existing CBCT rigid motion estimation methods rely on non-gradient optimization, which not only easily falls into local minimum, but also is very time-consuming.

[0005] Therefore, in view of the deficiencies of the prior art, it is necessary to provide a CBCT rigid motion artifact correction method and system based on a 3D-2D differentiable registration framework to solve the deficiencies of the prior art. SUMMARY

[0006] The first object of the present application is to provide a CBCT rigid motion artifact correction method based on a 3D-2D differentiable registration framework to avoid the deficiencies of the prior art. The CBCT rigid motion artifact correction method estimates the rigid disturbance that occurs during CBCT scanning, including mechanical motion, human voluntary / involuntary motion, by registering the estimated reconstructed image to the measured projection, and uses the estimated rigid disturbance to adjust the imaging geometry to reconstruct more accurate images, which are used to achieve more accurate rigid disturbance estimation.

[0007] The above object of the present application is achieved by the following technical measures:

[0008] A CBCT rigid motion artifact correction method based on a 3D-2D differentiable registration framework is provided, which is performed by the following steps:

[0009] S1, acquire CT images x with rigid motion artifacts, corresponding measurement projections y measure , CT geometry A geo and gantry rotation angles, the CT images x have a plurality of CT scan angles, and all the CT scan angles are arranged in order from front to back, and the learnable rigid motion parameters of the CT images x at each CT scan angle are defined, and the rigid motion parameters at all CT scan angles are initialized without motion;

[0010] S2, according to the CT geometry A geo and gantry rotation angles of S1, construct a forward projection sampling point grid at each CT scan angle, the forward projection sampling point grid is a set of sampling points, each sampling point is represented by (x, y, z) coordinates;

[0011] S3, at each CT scan angle, map the rigid motion parameters in S1 to a rigid transformation matrix by a rigid motion modeling method;

[0012] S4, at each CT scan angle, use the rigid transformation matrix of S3 to perform rigid spatial position transformation on the forward projection sampling grid of S2, to obtain a transformed forward projection sampling grid;

[0013] S5, at each CT scan angle, use the transformed forward projection sampling grid of S4 to perform differentiable interpolation sampling on the CT images x of S1, and perform sampling interval weighting and along-ray direction accumulation, to obtain a re-projection;

[0014] S6, at each CT scan angle, according to the re-projection of S5 and the measurement projection y measure of S1, calculate a cost function, and then update the rigid motion parameters by gradient backpropagation;

[0015] S7, define the first arranged CT scan angle as a reference angle, calculate the relative rigid disturbance of all CT scan angles relative to the reference angle, to obtain a plurality of relative rigid disturbances, and then update the rigid motion parameters at each CT scan angle using all the relative rigid disturbances;

[0016] S8, according to all the relative rigid disturbances obtained by S7, combine the CT geometry A geo of S1 to update, to obtain an updated CT geometry A geo ; and then reconstruct the CT image by a reconstruction algorithm and assign it to the CT image x, to obtain an updated CT image x;

[0017] S9, judge whether the iteration condition is reached, when yes, output the current updated CT image x and the rigid disturbance amount of all CT scan angles; when no, return to S3.

[0018] In the S1, the rigid motion parameters are defined as 3 translation amplitudes (t x , t y , t z ) and 3 Euler rotation angles (r x , r y , r z ) along three coordinate axes (x, y, z) in the world coordinate system.

[0019] In the S1, the rigid motion parameters are initialized without motion by setting t x , t y , t z , r x , r y and r z to 0 at each CT scan angle.

[0020] Preferably, the S2 comprises the following steps:

[0021] S2.1, according to the CT geometry A geo and the gantry rotation angle a j at the jth CT scan angle, obtaining the (x, y, z) coordinates of the X-ray source at the jth CT scan angle as S j , the (x, y, z) coordinates of the detector center as D j , and obtaining the detector row direction vector and the column direction vector

[0022] S2.2, according to the detector with the size (d row , d col ), obtaining the spatial position of each detector unit at the jth CT scan angle, which is represented by formula (1):

[0023]

[0024] wherein i is the ith row of the detector, i ranges from 1 to N, N is the number of rows of the detector, k is the kth column of the detector, k ranges from 1 to M, M is the number of columns of the detector, p det (i, k) is the detector unit with the row and column index (i, k);

[0025] S23, setting the sampling method as equal interval sampling and sampling L points on each line connecting the detector unit and the X-ray source, and constructing all the sampling points as the sampling G j at the jth CT scan angle, and G j ∈ R N×M×L×4 .

[0026] In the S3, the rigid transformation matrix is represented by equation (2):

[0027]

[0028] where T is a rigid transformation matrix, where R x is a rotation transformation matrix around the x-axis, R y is a rotation transformation matrix around the y-axis, R z is a rotation transformation matrix around the z-axis, and t is a translation vector.

[0029] Preferably, the rotation transformation matrix R x , the rotation transformation matrix R y , the rotation transformation matrix R z , and the translation vector t are represented by equations (3)-(5):

[0030]

[0031] Preferably, the S5 comprises the following steps:

[0032] S5.1, at each CT scan angle, using the forward projection sampling grid after transformation in S4 to perform differentiable interpolation sampling on the CT image x in S1;

[0033] S5.2, using trilinear interpolation to obtain the attenuation coefficient tensor μ j ∈ R N×M×L corresponding to the sampling grid after transformation at the jth CT scan angle;

[0034] S5.3, accumulating along the last sampling dimension of the attenuation coefficient tensor and multiplying by the sampling interval to obtain the back projection at the jth CT scan angle, which is represented by equation (7):

[0035]

[0036] where, is the signal value of the detector element p det (i, k) at the jth CT scan angle, is the sampling interval weight of the line connecting the detector element with the X-ray source at the row-column index (i, k), is the signal value of the detector element p det (i, k) at the jth CT scan angle, j is the attenuation value at the lth sampling point position on the line connecting the detector element with the X-ray source S

[0037] The cost function in the S6 is an L1 norm or an L2 norm.

[0038] The gradient rotation in the S6 is an implementation based on a pytorch framework.

[0039] Preferably, the reconstruction algorithm is an FDK algorithm, an algebraic reconstruction algorithm or a statistical reconstruction algorithm.

[0040] In the S9, the iteration condition is a maximum iteration number or a cost function value less than a threshold.

[0041] Preferably, the CT geometry A geo It comprises a parameter describing the spatial position relationship of the X-ray source, the scanned object and the detector.

[0042] Preferably, the CT geometry A geo It is a source-detector distance DSD, a source-rotation center distance, a detector unit size, a detector tilt angle and a detector array number.

[0043] The second object of the application is to provide a CBCT rigid motion artifact correction system based on a 3D-2D differentiable registration framework to avoid the defects of the prior art. The CBCT rigid motion artifact correction system based on the 3D-2D differentiable registration framework estimates the rigid disturbance amount, including mechanical motion, human voluntary / involuntary motion, occurring in the CBCT scanning process by registering the estimated reconstructed image to the measured projection, and the estimated rigid disturbance amount is used to adjust the imaging geometry to reconstruct more accurate images, and these images are used to realize more accurate rigid disturbance estimation.

[0044] The above object of the application is achieved by the following technical measures:

[0045] The application provides a CBCT rigid motion artifact correction system based on a 3D-2D differentiable registration framework to perform the CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework.

[0046] The CBCT rigid motion artifact correction method and system based on the 3D-2D differentiable registration framework comprises the following steps: S1, acquiring a CT image x with rigid motion artifacts, corresponding measured projection y measure , CT geometry A geo and gantry rotation angle, the CT image x has a plurality of CT scan angles, and all the CT scan angles are arranged in order from front to back, and the learnable rigid motion parameters of the CT image x at each CT scan angle are defined, and the rigid motion parameters at all CT scan angles are motionless initialized; S2, according to the CT geometry Ageo and the gantry rotation angle, construct a forward projection sampling point grid at each CTCT scan angle, the forward projection sampling point grid being a set of sampling points, each sampling point being represented by (x, y, z) coordinates;

[0047] S3, at each CT scan angle, map the rigid motion parameters in S1 into a rigid transformation matrix by a rigid motion modeling method; S4, at each CT scan angle, perform a rigid spatial position transformation on the forward projection sampling grid in S2 using the rigid transformation matrix in S3 to obtain a transformed forward projection sampling grid; S5, at each CT scan angle, perform a differentiable interpolation sampling on the CT image x in S1 using the transformed forward projection sampling grid in S4, and perform a sampling interval weighting and a ray direction accumulation to obtain a re-projection; S6, at each CT scan angle, according to the re-projection in S5 and the measurement projection y in S1 measure , calculate a cost function, then update the rigid motion parameters by gradient backpropagation; S7, define the CT scan angle ranked first as a reference angle, calculate the relative rigid disturbance of all CT scan angles relative to the reference angle to obtain a plurality of relative rigid disturbances, then update the rigid motion parameters at each CT scan angle using all the relative rigid disturbances; S8, according to all the relative rigid disturbances obtained in S7, combine the CT geometry A geo in S1 to update the CT geometry A geo ; then reconstruct the CT image by a reconstruction algorithm and assign it to the CT image x to obtain an updated CT image x; S9, judge whether the iteration condition is reached, when yes, output the current updated CT image x and the rigid disturbance amount of all CT scan angles; when no, return to S3. The present application estimates the rigid disturbance amount occurring in the CBCT scanning process, including mechanical motion, human voluntary / involuntary motion, by registering the estimated reconstructed image to the measurement projection, and the estimated rigid disturbance amount is used to adjust the imaging geometry to reconstruct a more accurate image, and these images are used to achieve more accurate rigid disturbance estimation. The beneficial effects of the present application are: 1, the method of the present application uses a rigid motion modeling method to strictly constrain the transformation matrix to satisfy the definition of rigid transformation. 2, the present application constructs a 3D-2D registration framework based on gradient, and estimates the rigid disturbance by registering the reconstructed image to the measurement projection. 3, the present application constrains the rigid disturbance to ensure that the image content is fixed in a specific imaging area, thereby ensuring the convergence of the rigid disturbance estimation. BRIEF DESCRIPTION OF DRAWINGS

[0048] The present application is further illustrated by the accompanying drawings, but the content in the drawings does not constitute any limitation on the present application.

[0049] Figure 1 It is a flowchart of a CBCT rigid motion artifact correction method based on a 3D-2D differentiable registration framework.

[0050] Figure 2 The result comparison chart before and after the algorithm correction of the application. DETAILED DESCRIPTION

[0051] The technical solutions of the application are further illustrated in combination with the following examples.

[0052] Example 1

[0053] A CBCT rigid motion artifact correction method based on a 3D-2D differentiable registration framework, such as Figure 1 , is performed by the following steps:

[0054] S1, obtaining a CT image x with rigid motion artifacts, corresponding measurement projection y measure , CT geometry A geo , and gantry rotation angle, while defining the learnable rigid motion parameters of the CT image x at each CT scan angle, and initializing the rigid motion parameters without motion;

[0055] S2, constructing a forward projection sampling point grid at each CT scan angle according to the CT geometry A geo and gantry rotation angle of S1, the forward projection sampling point grid being a set of sampling points, each sampling point being represented by (x, y, z) coordinates;

[0056] S3, at each CT scan angle, mapping the rigid motion parameters in S1 to a rigid transformation matrix through a rigid motion modeling method;

[0057] S4, at each CT scan angle, using the rigid transformation matrix of S3 to perform a rigid spatial position transformation on the forward projection sampling grid of S2, to obtain a transformed forward projection sampling grid;

[0058] S5, at each CT scan angle, using the transformed forward projection sampling grid of S4 to perform differentiable interpolation sampling on the CT image x of S1, and performing sampling interval weighting and along-ray direction accumulation, to obtain a re-projection;

[0059] S6, at each CT scan angle, calculating a cost function according to the re-projection of S5 and the measurement projection y measure of S1, and then updating the rigid motion parameters through gradient backpropagation;

[0060] S7, defining the CT scan angle with the first arrangement as the reference angle, calculating the relative rigid disturbance of all CT scan angles relative to the reference angle, to obtain a plurality of relative rigid disturbances, and then updating the rigid motion parameters at each CT scan angle using all the relative rigid disturbances;

[0061] S8, combine all the relative rigid perturbations obtained in S7 with the CT geometry A of S1 geo Update to obtain the updated CT geometry A geo ; and then reconstruct the CT image by the reconstruction algorithm and assign it to the CT image x to obtain the updated CT image x.

[0062] S9, determine whether the iteration condition is reached, and output the current updated CT image x and the rigid perturbation amount of all CT scan angles when yes; and return to S3 when no.

[0063] It should be noted that the relative rigid perturbation in S8 combines the imaging geometry to update the CT geometry, which means that the rigid motion of the object scanning part is converted into the spatial pose transformation of the source-detector pair, including modifying the source-detector position and direction. The measured projection y of the present application measure refers to the true projection data of the CT scanning part of the object.

[0064] In S1, the rigid motion parameters are defined as 3 translation amplitudes (t x , t y , t z ) and 3 Euler rotation angles (r x , r y , r z ) along the three coordinate axes (x, y, z) in the world coordinate system; in S1, the rigid motion parameters are initialized without motion, that is, t x , t y , t z , r x , r y and r z are all set to 0 at each CT scan angle.

[0065] It should be noted that the rigid motion parameters of the present application are used to describe the rigid perturbation of the object scanning part, and the specific form of the rigid motion parameters is determined according to the rigid motion modeling method adopted, including but not limited to the use of 3 translation amplitudes and 3 Euler rotation angles to describe the rigid motion. The motionless initialization of the present application means that the initialized rigid motion parameters describe that the object does not move during the scanning process.

[0066] Wherein, S2 includes the following steps:

[0067] S2.1, according to the CT geometry A geo of S1 and the gantry rotation angle α j of the jth CT scan angle, obtain the (x, y, z) coordinates of the X-ray source at the jth CT scan angle as S j , the (x, y, z) coordinates of the detector center as D j , and obtain the detector row direction vector and column direction vectors

[0068] S2.2, a detector with a detector unit size of (d row , d col ), the spatial position of each detector unit at the jth CT scan angle is obtained, and the spatial position is represented by formula (1):

[0069]

[0070] Wherein, i is the ith row of the detector, and i takes a value in the range of 1 to N, N is the number of rows of the detector, k is the kth column of the detector, k takes a value in the range of 1 to M, M is the number of columns of the detector, p det (i, k) is the detector unit of the row and column index (i, k);

[0071] S23, set the sampling mode to equal interval sampling and sample L points on each connecting line between the detector unit and the X-ray source, and all the sampling points are constructed as the sampling G j , and G j ∈R N×M×L×4 In order to facilitate spatial transformation, the spatial coordinates of the sampling points are represented in homogeneous form.

[0072] In S3, the rigid transformation matrix is represented by formula (2):

[0073]

[0074] Wherein, T is a rigid transformation matrix, wherein R x is a rotation transformation matrix around the x-axis, R y is a rotation transformation matrix around the y-axis, R z is a rotation transformation matrix around the z-axis, and t is a translation vector;

[0075] The rotation transformation matrix R x , the rotation transformation matrix R y , the rotation transformation matrix R z , and the translation vector t are represented by formula (3)-formula (5):

[0076]

[0077] It should be noted that the sampling mode of the present application includes but is not limited to equal interval sampling.

[0078] Wherein, S5 includes the following steps:

[0079] S5.1, at each CT scan angle, the CT image x of S1 is differentially interpolated and sampled using the front projection sampling grid transformed by S4;

[0080] S5.2, using trilinear interpolation, the attenuation coefficient tensor μ corresponding to the transformed sampling grid under the jth CT scan angle is obtained j ∈R N×M×L ;

[0081] S5.3, along the last sampling dimension of the attenuation coefficient tensor, accumulation and multiplication by the sampling interval are performed, and the re-projection under the jth CT scan angle is obtained, which is represented by formula (7):

[0082]

[0083] wherein, is the signal value of the detector unit p det (i,k) under the jth CT scan angle, is the sampling interval weight of the line connecting the detector unit with the X-ray source under the row and column index (i,k), is the signal value of the detector unit p det (i,k) under the jth CT scan angle, j is the attenuation value at the lth sampling point position on the line connecting the detector unit with the X-ray source S, and L is the total sampling points of each line connecting the detector unit with the X-ray source, and l is the sampling point index.

[0084] It should be noted that in formula (7), that is, the sampling interval is weighted, is accumulated along the ray direction.

[0085] The cost function in S6 is L1 norm and L2 norm; the gradient rotation in S6 is an implementation based on the pytorch framework. The reconstruction algorithm is FDK algorithm, algebraic reconstruction algorithm or statistical reconstruction algorithm.

[0086] In S9, the iteration condition is the maximum number of iterations or the cost function value is less than the threshold. The maximum number of iterations is 500, that is, the iteration is stopped after 500 times.

[0087] CT geometry A geo includes parameters describing the spatial position relationship of the X-ray source, the scanned object and the detector. The CT geometry A geo is the source-detector distance DSD, the source-rotation center distance, the detector unit size, the detector tilt angle and the detector array number.

[0088] The CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework estimates the rigid disturbance amount, including mechanical motion, human voluntary / involuntary motion, that occurs in the CBCT scanning process by registering the estimated reconstructed image to the measured projection, and the estimated rigid disturbance amount is used to adjust the imaging geometry to reconstruct more accurate images, and these images are used to achieve more accurate rigid disturbance estimation. The beneficial effects of the present application are: 1. The method of the present application uses a rigid motion modeling method, which strictly constrains the transformation matrix to satisfy the definition of rigid transformation. 2. The present application constructs a 3D-2D registration framework based on gradient, and estimates the rigid disturbance by registering the reconstructed image to the measured projection. 3. By constraining the rigid disturbance, the present application ensures that the image content is fixed in a specific imaging area, thereby ensuring the convergence of the rigid disturbance estimation.

[0089] Embodiment 2

[0090] A CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework of embodiment 1 is applied. In this embodiment,

[0091] The CT image x0 of the patient's head scanned by the Fussen commercial dental CBCT (Fussen Matrix5000) is used to simulate the case of rigid motion of the patient's head during CT scanning. Specifically, the motion amplitude values of t x , t y , t z , r x , r y and r z are randomly generated at each CT scanning angle, combined with the CT geometry A geo , the corresponding forward projection geometry is generated, and the CT image x0 is forward projected to simulate the measured projection y measure obtained under the condition of rigid motion of the patient's head during CT scanning.

[0092] Without considering the motion condition, the CT image x with motion artifacts is directly reconstructed using the CT geometry A geo . The CT image x with motion artifacts, the measured projection y measure and the CT geometry A geo are input into the CBCT rigid motion artifact correction algorithm based on the 3D-2D differentiable registration framework, and the algorithm is run to estimate the rigid motion parameters and compensate for the motion, and finally the corrected CT image and the rigid disturbance amount at all CT scanning angles are obtained.

[0093] Specifically, the CT geometry A geoThe image size of the CT image is set to 800x800x400, the source-to-detector distance (DSD) is set to 635 mm, the source-to-rotation center distance (DSO) is set to 425 mm, the size of the detector cell is 0.2 mmx0.2 mm, the detector array is 500x1500, the reconstructed image pixel size is 0.2 mmx0.2 mmx0.2 mm, and the number of projections is 300.

[0094] In S6, the re-projection of S5 and the measured projection y of S1 are used to calculate the rigid motion parameters at each CT scan angle. measure The cost function is calculated, and then the gradient backpropagation is used to update the rigid motion parameters. In this embodiment, the cost function is specifically an L1 norm of the re-projection and the measured projection, and the cost function Loss is represented by the following formula:

[0095]

[0096] wherein, is the re-projection of the detector cell with the row and column indexes (i, k) at the jth CT scan angle, is the measured projection y of the detector cell with the row and column indexes (i, k) at the jth CT scan angle, measure J is the number of projection angles, N is the number of rows of the detector, and M is the number of columns of the detector.

[0097] The gradient backpropagation is realized based on the pytorch framework. It should be noted that the gradient backpropagation is a common operation in the technical field, and will not be described here.

[0098] In S9, the iteration criterion is that the maximum number of iterations is reached, and the maximum number of iterations is set to 500.

[0099] In this example, the CBCT rigid motion artifact correction effect based on the 3D-2D differentiable registration framework is demonstrated. Figure 2 The CBCT rigid motion artifact correction effect based on the 3D-2D differentiable registration framework is demonstrated.

[0100] The CBCT rigid motion artifact correction effect based on the 3D-2D differentiable registration framework is demonstrated. Figure 2 It can be seen that the rigid motion artifact CT image has obvious double-contour artifacts and distorted anatomical structures. The present application can effectively restore details such as teeth and bone tissue from the severely distorted rigid motion artifact CT image and eliminate the double-contour artifacts. This result verifies the feasibility and effectiveness of the CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework proposed in the present application.

[0101] Embodiment 3

[0102] A system for CBCT rigid motion artifact correction based on a 3D-2D differentiable registration framework, wherein the CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework of embodiment 1 is performed.

[0103] The system of CBCT rigid motion artifact correction based on the 3D-2D differentiable registration framework estimates the rigid disturbance amount, including mechanical motion, human voluntary / involuntary motion, occurring in the CBCT scanning process by registering the estimated reconstructed image to the measured projection, and the estimated rigid disturbance amount is used to adjust the imaging geometry to reconstruct more accurate images, and these images are used to achieve more accurate rigid disturbance estimation. The beneficial effects of the present application are: 1. The method of the present application adopts a rigid motion modeling method, which strictly constrains the transformation matrix to satisfy the definition of rigid transformation. 2. The present application constructs a 3D-2D registration framework based on gradient, and estimates the rigid disturbance by registering the reconstructed image to the measured projection. 3. The present application ensures the convergence of the rigid disturbance estimation by constraining the rigid disturbance and ensuring that the image content is fixed in a specific imaging area.

[0104] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application and not to limit the protection scope of the present application. Although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present application can be modified or replaced by equivalents without departing from the essence and scope of the technical solutions of the present application.

Claims

1. A CBCT rigid motion artifact correction method based on a 3D-2D differentiable registration framework, characterized in that, The method comprises the following steps: S1, acquire a CT image x with a rigid motion artifact, a corresponding measurement projection y measure , a CT geometry A geo , and a gantry rotation angle, the CT image x has a plurality of CT scan angles, and all the CT scan angles are arranged in order from front to back, and the learnable rigid motion parameters of the CT image x at each CT scan angle are defined, and the rigid motion parameters at all CT scan angles are initialized without motion; S2. The CT geometry A according to S1 geo and the gantry rotation angle, construct a forward projection sample point grid at each CT scan angle, the forward projection sample point grid being a set of sample points, each sample point being represented by a (x, y, z) coordinate; S3. mapping the rigid motion parameters in S1 to a rigid transformation matrix by a rigid motion modeling method at each CT scan angle; S4. performing a rigid spatial position transformation on the forward projection sampling grid in S2 using the rigid transformation matrix in S3 at each CT scan angle to obtain a transformed forward projection sampling grid; S5. performing a differentiable interpolation sampling on the CT image x in S1 using the transformed forward projection sampling grid in S4 at each CT scan angle, and performing a sampling interval weighting and a ray direction accumulation to obtain a re-projection; S6. At each CT scan angle, the re-projection according to S5 and the measured projection y according to S1 measure The cost function is calculated, and then the rigid motion parameters are updated by gradient backpropagation. S7. defining the CT scan angle with the first rank as a reference angle, calculating relative rigid perturbations of all CT scan angles relative to the reference angle to obtain a plurality of relative rigid perturbations, and then updating the rigid motion parameters at each CT scan angle using all the relative rigid perturbations; S8. All relative rigid perturbations resulting from S7 combined with the CT geometry A of S1 geo An update is performed to obtain an updated CT geometry A geo ; and the CT image is reconstructed by a reconstruction algorithm and assigned to the CT image x to obtain an updated CT image x; S9. determining whether an iteration condition is reached, outputting the current updated CT image x and the rigid perturbation of all CT scan angles when the iteration condition is reached, and returning to S3 when the iteration condition is not reached.

2. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of claim 1, wherein: In the S1, the rigid motion parameters are defined as 3 translation amplitudes (t x , t y , t z ) and 3 Euler rotation angles (r x , r y , r z ) along the three coordinate axes (x, y, z) in the world coordinate system; In the S1, the rigid motion parameters are initialized to zero t x , t y , t z , r x , r y and r z are all set to 0.

3. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of claim 1, wherein, The S2 comprises the following steps: S2.

1. CT geometry A according to S1 geo and the gantry rotation angle a of the jth CT scan angle j , the (x, y, z) coordinates of the X-ray source at the jth CT scan angle is S j , the (x, y, z) coordinates of the detector center is D j , and the detector row direction vector is and the column direction vector is S2.2, according to the detector unit size (d row , d col ) of the detector, the spatial position of each detector unit at the jth CT scan angle is obtained, and the spatial position is represented by formula (1) where i is the i-th row of the detector, i ranges from 1 to N, N is the number of rows of the detector, k is the k-th column of the detector, k ranges from 1 to M, M is the number of columns of the detector, p det (i,k) is a detector cell at row-column index (i,k). S23, set the sampling mode as equal interval sampling and sample L points on each connecting line, sample each connecting line between the detector unit and the X-ray source, and all sampling points are constructed as the sampling G of the jth CT scanning angle j , and G j ∈R N×M×L×4 .

4. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of claim 1, wherein: In the S3, the rigid transformation matrix is represented by formula (2): where T is a rigid transformation matrix, where R x is a rotation transformation matrix about the x-axis, R y is a rotation transformation matrix about the y-axis, R z is a rotation transformation matrix about the z-axis, and t is a translation vector. Rotation transformation matrix R x Rotation transformation matrix R y Rotation transformation matrix R z Translation vector t is represented by equations (3) - (5):

5. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of claim 1, wherein, The S5 comprises the following steps: S5.

1. performing a differentiable interpolation sampling on the CT image x in S1 using the transformed forward projection sampling grid in S4 at each CT scan angle; S5.2, using trilinear interpolation, get the attenuation coefficient tensor μ corresponding to the sampling grid after transformation under the jth CT scan angle j ∈R N×M×L ; S5.

3. performing an accumulation along the last sampling dimension of the attenuation coefficient tensor and multiplying by the sampling interval to obtain a re-projection under the jth CT scan angle, the re-projection being represented by formula (7): wherein is the signal value of the detector cell p at row index (i, k) for the j-th CT scan angle det is the signal value of the detector cell p at row index (i, k) for the j-th CT scan angle is the sampling interval weight of the line between the detector cell at row index (i, k) and the X-ray source Sjfor the j-th CT scan angle is the signal value of the detector cell p at row index (i, k) for the j-th CT scan angle det is the signal value of the detector cell p at row index (i, k) for the j-th CT scan angle j is the attenuation value at the l-th sampling point position on the line between the detector cell at row index (i, k) and the X-ray source Sjfor the j-th CT scan angle, L is the total number of sampling points of each line between a detector cell and the X-ray source, and l is the sampling point index.

6. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of any one of claims 1 to 5, characterized in that: The cost function in the S6 is an L1 norm or an L2 norm; The gradient rotation in the S6 is an implementation based on a pytorch framework; The reconstruction algorithm is an FDK algorithm, an algebraic reconstruction algorithm or a statistical reconstruction algorithm.

7. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of any one of claims 1 to 5, characterized in that: In the S9, the iteration condition is a maximum iteration number or a cost function value less than a threshold.

8. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of any one of claims 1 to 5, characterized in that: The CT geometry A geo includes parameters that describe the spatial relationship of the positions of the x-ray source, the scanned object, and the detector.

9. The CBCT rigid motion artifact correction method based on 3D-2D differentiable registration framework of any one of claims 1 to 5, characterized in that: The CT geometry A geo are the source-detector distance DSD, the source-rotation center distance, the detector cell size, the detector tilt angle, the number of detector arrays.

10. A system for CBCT rigid motion artifact correction based on a 3D-2D differentiable registration framework, characterized in that: The CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework is executed. The CBCT rigid motion artifact correction method based on the 3D-2D differentiable registration framework is executed.

Images