Double-bulb tube ray source cone beam CT geometric correction and registration method based on spiral steel ball die body

Abstract

The invention discloses a double-bulb tube ray source cone beam CT (Computed Tomography) geometric correction and registration method, which specifically comprises the following steps of: carrying out double-ray source CT data acquisition on a spiral steel ball die body, and determining an approximate light source vertical incidence center position according to a maximum-density projection image of a projection drawing; estimating ideal geometry and obtaining a three-dimensional CT image of the spiral steel ball die body corresponding to the dual-ray source through a filtering back projection algorithm; recognizing and recording steel ball positions in the two three-dimensional CT images; matching the steel balls in the two three-dimensional CT images according to the theoretical height difference of the double-ray source, and recording the position containing the steel balls together; fitting the real image coordinate groups of the two sources to the same coordinate system through rigid body transformation according to the steel balls contained in the two sources; and substituting projection matrix reconstruction to carry out filtering back projection reconstruction on the dual-ray source CT image. Through the steps, the geometric artifacts of the double-source CT image can be effectively removed, image registration of the double-source CT can be completed, and the high-quality double-source reconstructed image after geometric registration can be output on real double-source CT data.

Description

Technical Field

[0001] The invention belongs to the technical field of CT image artifact correction and registration, and in particular relates to a double-tube ray source cone-beam CT geometric correction and registration method based on a spiral steel ball phantom. Background Art

[0002] In medical imaging, cone-beam CT equipped with dual-tube radiation sources arranged along the Z direction can significantly increase the imaging field of view in the Z direction and reduce cone-beam artifacts. Due to machine jitter during the scanning process, CT images will produce geometric artifacts, which appear radially on the image and will significantly affect the diagnostic accuracy and quality of the image. The images reconstructed by two X-ray sources need to be calibrated to ensure that the three-dimensional CT images are accurately aligned in the same coordinate system to ensure the clarity and accuracy of the stitched image. Summary of the invention

[0003] The technical problem to be solved by the present invention is to provide a dual-source CT geometric artifact removal and geometric registration method based on offline correction in order to solve the above-mentioned problem of removing geometric artifacts and completing geometric registration at the same time.

[0004] In order to achieve the above technical objectives, the technical solutions adopted by the present invention are as follows:

[0005] A dual-source CT geometric artifact removal and geometric registration method based on offline correction, the method comprising the following steps:

[0006] Step 1: Perform dual-ray source CT data acquisition on the spiral steel ball phantom, determine the approximate vertical incidence center position of the light source based on the maximum density projection image, estimate the ideal geometry, and obtain the three-dimensional CT image of the spiral steel ball phantom corresponding to the dual-ray source through the filtered back projection algorithm;

[0007] Step 2: Identify and record the positions of the steel balls in the two 3D CT images to obtain two sets of steel ball coordinates;

[0008] Step 3: According to the theoretical height difference of the dual-ray source, match the steel balls in the two 3D CT images and record the positions of the steel balls;

[0009] Step 4: Fit the steel ball coordinate groups of the two sources to the same coordinate system through rigid body transformation based on the common steel balls;

[0010] Step 5: Match the three-dimensional coordinates of the steel ball and the two-dimensional coordinates on the corresponding projection map for the dual-ray sources through ideal geometric projection;

[0011] Step 6: Combine the three-dimensional coordinates and projection of the steel ball in the common coordinate system Figure 2 dimensional coordinates to solve the projection matrix;

[0012] Step 7: Substitute the projection matrix into the reconstruction to perform filtered back-projection reconstruction on the dual-ray source CT image.

[0013] To optimize the above technical solutions, the specific measures taken also include:

[0014] In step 1, the spiral steel ball model is Dual-ray source After CT data acquisition, on the projection image, the center of the vertical incident center of the light source is in the middle of the row, but the center of the column needs to pass Determination of the trajectory of the steel ball For the maximum density projection image, the motion trajectory of the steel ball is obtained by taking the maximum value of the z-axis. The column center position of the detector is determined by the motion trajectory of the steel ball.

[0015] After subtracting half of the projected image column pixel from the determined column center, multiply it by the physical size of the pixel value to obtain the detector offset in this direction. Perform the above operation on both sources to obtain the approximate vertical incidence center of the two ray sources, which are recorded as S1 and S2.

[0016] In step 2, for identifying the steel balls in the reconstructed image, after threshold processing is performed on the reconstructed image, each steel ball appears as a connected component in the binary image, and the center position of each connected component is recorded. The steel ball positions of the images reconstructed by light sources #1 and #2 are grouped and recorded as K1 and K2 respectively.

[0017] In step 5, the three-dimensional coordinates of the steel ball need to be projected using the projection matrix of ideal geometry. Substituting the 3-column and 4-row projection matrix of ideal geometry can obtain the two-dimensional projection coordinates (u, v) of the detector plane corresponding to the three-dimensional spatial coordinates (x, y, z) of each steel ball.

[0018] Compared with the prior art, the advantages of the present invention are as follows: the present invention can effectively remove geometric artifacts of dual-source CT images and complete image registration of dual-source CT, and can output high-quality dual-source reconstructed images after geometric registration on real dual-source CT data. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 It is a schematic diagram of the process of the present invention;

[0020] Figure 2 This is an example diagram of the motion trajectory of the steel ball obtained in step 1;

[0021] Figure 3 is S1, S2, S in step 4 mid Example diagram of the relationship of three points in the detector plane;

[0022] Figure 4 This is an example of the curve drawing results of T1 and T2 in step 4;

[0023] Figure 5This is an example diagram of calculating the projection matrix in three-dimensional space;

[0024] Figure 6 These are examples of pig feet projection images reconstructed using projection matrices (3D perspective);

[0025] Figure 7 This is an example of the output result after the pig's foot data is geometrically registered (3D perspective);

[0026] Figure 8 It is the correspondence diagram between the three-dimensional coordinates of the steel ball and the two-dimensional coordinates on the projection diagram. DETAILED DESCRIPTION

[0027] The embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings.

[0028] Example: Figure 1 As shown, the present invention is a double-tube ray source cone beam CT geometric correction and registration method based on a spiral steel ball phantom, comprising the following steps:

[0029] Step 1: Number the steel ball points in the spiral model as P n , dual-ray source CT data acquisition was performed on the spiral phantom, and the approximate vertical incidence center of the light source was determined based on the maximum density projection image. The ideal geometry was estimated and the three-dimensional CT image of the spiral steel ball phantom corresponding to the dual-ray source was obtained through the filtered back-projection algorithm.

[0030] Specifically: the spiral steel ball phantom Dual-ray source After CT data acquisition, on the projection image, the center of the vertical incident center of the light source is in the middle of the row, but the center of the column needs to pass Determination of the trajectory of the steel ball For the maximum density projection image, the motion trajectory of the steel ball is obtained by taking the maximum value of the z-axis. The column center position of the detector is determined by the motion trajectory of the steel ball.

[0031] After subtracting half of the projected image column pixel from the determined column center, multiply it by the physical size of the pixel value to obtain the detector offset in this direction. Perform the above operation on both sources to obtain the approximate z-direction coordinates of the vertical incident center of the light source of the two sources, recorded as z1 and z2.

[0032] Step 2: Identify the steel balls from the two sets of reconstruction images and obtain the coordinate sets K1 and K2 of the steel balls; observe the steel ball number P of the spiral phantom n The corresponding relationship between the steel ball points in the steel ball coordinate group is recorded as P, and the steel ball number actually represented by each coordinate in the steel ball coordinate group K1, K2 is recorded as P 1,n , P 2,n , recorded in P 1,n and P 2,n The steel balls that have appeared are numbered as P 3,n .

[0033] The specific steps of identifying the steel balls in the reconstructed image are as follows: threshold processing is performed on the original 3D reconstructed image. By selecting a suitable threshold, the steel ball area in the image is extracted to generate a 3D binary image. The function is used to identify all connected areas in the image, and the center position of each connected area is obtained to obtain the center coordinates of each steel ball.

[0034] Step 3: According to the theoretical height difference of the dual-ray source, match the steel balls in the two 3D CT images and record the positions of the steel balls;

[0035] Specifically, for a steel ball in the phantom, its x and y coordinates are basically consistent in the two ray source reconstruction images. In the z direction, the difference in its coordinates (physical size) is the height difference L = |S1-S2| between S1 and S2.

[0036] After adding L to the z coordinate of the steel ball in K2, sort K1 and K2 by z coordinate, and find the corresponding coordinate in K1 from the first steel ball in K2 through nearest neighbor matching. Because they have been sorted by z coordinate, the matching is sequential in K1 and K2. When the matching reaches the end of K1, the matching ends and the common steel balls are obtained. The coordinates of the common steel balls recorded in K1 and K2 are Rc1 and Rc2 respectively.

[0037] Step 4: The coordinate system centers of the steel ball coordinate groups K1 and K2 are respectively located at the approximate vertical incidence centers of the two ray sources z1 and z2. The steel ball coordinate groups of the two sources are fitted to the same coordinate system through rigid body transformation based on the common steel balls, and the coordinate system center is the midpoint S1 and S2. mid . S1, S2, S mid The relationship between the three points on the detector plane is as follows: Figure 3 shown.

[0038] Specifically, take the average value of each corresponding steel ball point pair of Rc1 and Rc2 to obtain the coordinate group Rc. The rigid body transformation matrix from Rc1 to Rc is recorded as M1, and the rigid body transformation matrix from Rc2 to Rc is recorded as M2. Perform a rigid body transformation on K1 to transfer its steel ball coordinates to the coordinate system of Rc. The result of M1*K1 is recorded as T1. Similarly, for K2, the result of M2*K2 is recorded as T2. The curve drawing results of T1 and T2 are as follows: Figure 4 , T1 is the blue curve, T2 is the black curve, and the overlapping part of the two is the red curve.

[0039] Step 5: For each View of the projection image of the first ray source, identify the coordinates of the steel ball points after threshold processing. Use the projection matrix of ideal geometry to perform orthographic projection on the steel ball coordinate group K1, and use the points obtained by orthographic projection to match the nearest neighbor steel ball points identified on the projection image for each steel ball, and obtain the corresponding steel ball two-dimensional coordinate group W1 of the steel ball coordinate group K1 on the projection image at each angle; repeat the above steps for the projection image of the second ray source and the steel ball coordinate group K2 to obtain the corresponding steel ball two-dimensional coordinate group W2 of the steel ball coordinate group K2 on the projection image at each angle.

[0040] as follows Figure 8 For a certain example, the black dots are the two-dimensional coordinates of the steel ball obtained by orthographic projection at a certain projection angle, and the blue circles are the identified coordinates of the steel ball on the projection image. The two will be matched by nearest neighbors to obtain the correspondence between the three-dimensional coordinates of the steel ball and the two-dimensional coordinates on the projection image.

[0041] Step 6: Combine T1 and W1 to obtain the projection matrix Pmatrix1, and combine T2 and W2 to obtain the projection matrix Pmatrix2.

[0042] Specifically, for a steel ball, the ideal distribution three-dimensional coordinate set in T is (x i ,y i ,z i ), the two-dimensional projection coordinates of the detector plane in W are recorded as (u i ,v i ). The 3-row 4-column matrix P represents the mapping relationship between the 3D space coordinates and the projection plane coordinates. i For the scaling parameter, this equation can be established:

[0043]

[0044] These equations can be written explicitly in terms of the elements of P as follows:

[0045]

[0046] The variants are:

[0047]

[0048] Multiple groups (x i ,y i ,z i ) and (u i ,v i ) Yes, we can build equations to solve for the elements of P. First, set P 34 =1, and then use the least squares method to solve for P.

[0049] Step 7:

[0050] The mapping relationship represented by Pmatrix1 and Pmatrix2 will map the projection coordinates of the two sources to the same three-dimensional space, with the coordinate center being S mid Substituting the projection matrix into the reconstruction can move the coordinate system centers of both reconstruction images to S mid , that is, the geometric registration is completed. At the same time, the overlapped parts of the two reconstructed images can be combined with the information of the two reconstructed images by weighted averaging to output the final reconstruction result.

[0051] Specifically:

[0052] For the data to be corrected, Pmatrix1 and Pmatrix2 are used to reconstruct the two source projection images to obtain the reconstructed images Re1 and Re2, as shown in Figure 6 . Check whether each pixel value of Re1 and Re2 is within the reconstruction range (i.e., HU value > -1000), thereby generating two masks mask_s1_o and mask_s2_o, which are used to mark the valid areas in the image. Add the two masks together to obtain a merged mask mask to identify the valid areas in the two reconstructed images. For pixel positions with a value of 2 in the merged mask (i.e., both masks are 1 at this position), adjust their values ​​to 0.5, thereby averaging the values ​​of the two reconstructed images at these positions. By applying the corresponding mask to each reconstructed image and adding the results, a combined reconstructed image is generated, as shown in Figure 7 .

[0053] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A dual-source CT geometric registration and geometric artifact removal method based on offline correction, characterized in that: The method comprises the following steps: Step 1: Perform dual-ray source CT data acquisition on the spiral steel ball phantom, determine the approximate vertical incidence center position of the light source based on the maximum density projection image, estimate the ideal geometry, and obtain the three-dimensional CT image of the spiral steel ball phantom corresponding to the dual-ray source through the filtered back-projection algorithm; Step 2: Identify and record the positions of the steel balls in the two 3D CT images to obtain two sets of steel ball coordinates; Step 3: According to the theoretical height difference of the dual-ray source, match the steel balls in the two 3D CT images and record the positions of the steel balls; Step 4: Fit the coordinate groups of the two sources into the same coordinate system through rigid body transformation based on the common inclusion of the steel ball; Step 5: Match the three-dimensional coordinates of the steel ball and the two-dimensional coordinates on the corresponding projection map for the dual-ray sources through ideal geometric projection; Step 6: Combine the three-dimensional coordinates of the steel ball space and the two-dimensional coordinates of the projection image in the common coordinate system to solve the projection matrix; Step 7: Substitute the projection matrix into the reconstruction to perform filtered back-projection reconstruction on the dual-ray source CT image.

2. The method for dual-source CT geometric registration and geometric artifact removal based on offline correction according to claim 1, characterized in that: In step 1, the steel ball points in the spiral model are numbered P n , dual-ray source CT data acquisition was performed on the spiral phantom, the approximate vertical incidence center of the light source was determined based on the maximum density projection image, the ideal geometry was estimated, and the three-dimensional CT image of the spiral steel ball phantom corresponding to the dual-ray source was obtained through the filtered back-projection algorithm, Specifically, after the dual-ray source CT data of the spiral steel ball phantom is collected, on the projection image, the center of the vertical incident center of the light source is in the middle of the row, and the center of the column needs to be determined by the motion trajectory of the steel ball. For the maximum density projection image, the motion trajectory of the steel ball is obtained by taking the maximum value of the z-axis, and the column center position of the detector is determined by the motion trajectory of the steel ball. After subtracting half of the projected image column pixels from the determined column center, multiply it by the physical size of the pixel value to obtain the offset of the detector in this direction. The above operation is performed on both sources to obtain the approximate z-direction coordinates of the vertical incident center of the light sources of the two sources, which are recorded as z1 and z2.

3. The method for dual-source CT geometric registration and geometric artifact removal based on offline correction according to claim 2, characterized in that: Step 2 is as follows: The steel balls are identified by two sets of reconstruction images and the coordinate sets K1 and K2 of the steel balls are obtained; the steel ball number P of the spiral phantom is observed. n The corresponding relationship between the steel ball points in the steel ball coordinate group is recorded as P, and the steel ball number actually represented by each coordinate in the steel ball coordinate group K1, K2 is recorded as P 1,n , P 2,n , recorded in P 1,n and P 2,n The steel balls that have appeared are numbered as P 3,n ,The specific steps for identifying the steel balls in the reconstructed image are: threshold processing the original 3D reconstructed image, extracting the steel ball area in the image by selecting a suitable threshold, generating a 3D binary image, using the function to identify all connected areas in the image, and obtaining the center position of each connected area to obtain the center coordinates of each steel ball.

4. The method for dual-source CT geometric registration and geometric artifact removal based on offline correction according to claim 3, characterized in that: The step 3 is as follows: For a steel ball in the phantom, the x and y coordinates are basically consistent when reconstructed by two ray sources. The biological coordinate difference in the z direction (physical size) is the height difference in the z direction of the approximate vertical incident centers of the two ray sources. Assuming that the z-direction coordinates corresponding to light sources #1 and #2 are z1 and z2 respectively, the coordinate difference L = |z1-z2|. After adding L to the z coordinate of the steel ball in the steel ball coordinate group K2 of the second ray source, sort the steel ball coordinate groups K1 and K2 of the two ray sources by z coordinate, and find the corresponding coordinate in K1 from the first steel ball coordinate of K2 through nearest neighbor matching. Because they have been sorted by z coordinate, the matching is sequential in K1 and K2. When the match reaches the end of K1, the matching ends and the commonly included steel balls are obtained. The coordinates of the commonly included steel balls recorded in K1 and K2 are Rc1 and Rc2 respectively.

5. The method for dual-source CT geometric registration and geometric artifact removal based on offline correction according to claim 4, characterized in that: The specific details of step 5 are as follows: for each projection viewing angle of the projection image of the first ray source, the coordinates of the steel ball points are identified after threshold processing, and the steel ball coordinate group K1 is orthographically projected using the projection matrix of ideal geometry. For each steel ball, the points obtained by the orthographic projection are used to match the nearest neighbor steel ball points identified on the projection image to obtain the steel ball two-dimensional coordinate group W1 corresponding to the steel ball coordinate group K1 on the projection image at each angle; repeat the above steps for the projection image of the second ray source and the steel ball coordinate group K2 to obtain the steel ball two-dimensional coordinate group W2 corresponding to the steel ball coordinate group K2 on the projection image at each angle.

6. The method for dual-source CT geometric registration and geometric artifact removal based on offline correction according to claim 5, characterized in that: Step 6: Combine T1 and W1 to obtain the projection matrix Pmatrix1, and combine T2 and W2 to obtain the projection matrix Pmatrix2. Specifically, for a steel ball, the ideal distribution three-dimensional coordinate set in T is (x i ,y i ,z i ), the two-dimensional projection coordinates of the detector plane in W are recorded as (u i ,v i ), the 3-row 4-column matrix P represents the mapping relationship between the 3D space coordinates and the projection plane coordinates, M i To scale the parameters, establish this equation: These equations can be written explicitly in terms of the elements of P as follows: The variants are: Obviously, we get multiple groups (x i ,y i ,z i ) and (u i ,v i ) Yes, we can set up equations to solve for the elements of P.

7. The method for dual-source CT geometric registration and geometric artifact removal based on offline correction according to claim 6, characterized in that: Step 7: The mapping relationship represented by Pmatrix1 and Pmatrix2 will map the projection coordinates of the two sources to the same three-dimensional space, with the coordinate center being S mid , substituting the projection matrix into the reconstruction can move the coordinate system centers of the two reconstructed images to S mid , that is, the geometric registration is completed. At the same time, the overlapping parts of the two reconstructed images can be combined with the information of the two reconstructed images by weighted averaging to output the final reconstruction result. Specifically: For the data to be corrected, Pmatrix1 and Pmatrix2 are used to reconstruct the two source projection images to obtain the reconstructed images Re1 and Re2 respectively. Each pixel value of Re1 and Re2 is checked to see if it is within the reconstruction range (i.e., HU value>-1000), thereby generating two masks mask_s1_o and mask_s2_o, which are used to mark the valid areas in the image. The two masks are added together to obtain a merged mask mask to identify the valid areas in the two reconstructed images. For the pixel positions with a value of 2 in the merged mask (i.e., both masks are 1 at this position), their values ​​are adjusted to 0.5, thereby averaging the values ​​of the two reconstructed images at these positions. By applying the corresponding mask to each reconstructed image and adding the results, a comprehensive reconstructed image is generated.

Images