Dispersive metal artifact removal method based on energy spectrum estimation technology

Abstract

A dispersive metal artifact removal method based on an energy spectrum estimation technology belongs to the field of CT imaging technology and artifact removal, and comprises the following steps: initializing an X-ray source target material and element excitation characteristic energy; a step-shaped phantom is placed; starting a ray source, and obtaining step phantom projection data; reversely estimating an X-ray source energy spectrum curve by using a conjugate gradient method; the phantom is removed, and the measured workpiece is scanned; and performing energy spectrum correction on a projected image obtained by scanning, performing reconstruction by using an FDK algorithm, and finally outputting a high-quality CT reconstruction image without metal artifacts. According to the method, the metal artifact phenomenon caused by beam hardening and photon starvation can be effectively reduced, operation is easy and convenient, redundant image post-processing algorithms are not needed, and the CT reconstruction quality can be remarkably improved.

Description

Technical Field

[0001] The present invention relates to the field of CT imaging and image processing, and more specifically, to a method for removing diffuse metal artifacts based on energy spectrum estimation technology. Background Art

[0002] In recent years, X-ray computed tomography has been widely used in the field of industrial nondestructive testing. By performing projection scans of the object at multiple different angles and reconstructing the projection data using computational methods, the three-dimensional internal structure information of the object can be obtained. However, when the object being inspected contains metal materials with high atomic numbers, beam hardening and photon starvation will cause severe light and dark striped metal artifacts, which seriously affect the image quality. In response to such problems, many patents have designed methods for removing metal artifacts from CT images based on threshold segmentation or deep learning. For example, the patent with publication number CN116543063A has designed a method for removing metal artifacts from CT images. By constructing an initial adaptive iterative learning model based on wavelet transform, the CT image is decomposed by using the target adaptive iterative learning model optimized by the optimization objective function, the area of ​​metal artifacts is calculated, and the metal artifacts in the CT image are removed to obtain a CT image with reduced metal artifacts. Another patent with publication number CN117934348A proposes a method and device for removing metal artifacts based on a generative adversarial network. By deeply training images with and without metal artifacts, and using a generator and a discriminator to build a model, this method effectively removes metal artifacts from the initial image and significantly improves the clarity and usability of the image after artifact removal. The patent with publication number CN118314230A discloses a method, system and device for removing metal artifacts. The method includes obtaining an original CT image, separating metal areas and filling corresponding areas, and finally obtaining a CT image with artifacts removed. Through this method, the impact of metal artifacts on image quality is significantly reduced, and the clarity and accuracy of CT images are improved, making medical imaging diagnosis more reliable. However, although the method of removing metal artifacts using post-image related technologies can eliminate metal artifacts to a certain extent and improve image quality, it is difficult to handle once it involves complex metal artifact types, such as diffuse metal artifacts with complex streaks. In addition, this type of method does not fundamentally obtain the real information obscured by metal artifacts, but uses estimation or speculation based on prior information, which lacks authenticity and cannot meet the requirements of industrial non-destructive testing for accuracy and quality.

[0003] Considering the essential reasons for beam hardening or photon starvation, and the fact that the actual X-ray energy spectrum is a spectrum of continuous distribution of multiple energies, while the reconstruction process of CT images idealizes it into a spectrum of single energy, the final image quality is poor. The patent with publication number CN116664429A designs a semi-supervised method for removing metal artifacts from multi-spectrum CT images. Based on the semi-supervised average teacher-student model, the adopted U-net network structure adds a slice feature attention module, increases the perception of deep feature information between channels, and improves the feature connection between different single-energy images of multi-spectrum CT slices, thus improving the image quality. Although this method introduces more information for the elimination of metal artifacts from the perspective of energy spectrum, it is still a metal artifact elimination method from the perspective of data image post-processing. Its process is complicated and requires additional data sets for training, which cannot meet the efficiency requirements of industrial non-destructive testing.

[0004] Therefore, we consider implementing a diffuse metal artifact removal method based on energy spectrum estimation technology. We use a phantom made of a material that matches the X-ray source tube voltage to estimate and measure the X-ray source energy spectrum, then perform a three-dimensional scan of the workpiece to be measured. At the same time, we use the approximate filling material and its proportion of the workpiece to correct the projection data. Finally, we output a high-quality CT reconstructed image without metal artifacts, which further improves the accuracy and efficiency of industrial nondestructive testing. Summary of the invention

[0005] The purpose of the present invention is to provide a diffuse metal artifact removal method based on energy spectrum estimation technology to solve the problem that metal artifacts in industrial CT images, especially diffuse metal artifacts, are difficult to remove. The present invention uses a phantom made of a material that matches the voltage of the X-ray source tube to estimate and measure the X-ray source energy spectrum, and then performs a three-dimensional scan on the measured workpiece. At the same time, the projection data is corrected using the approximate filling material and its proportion of the measured workpiece. The process is simple and efficient, and has strong feasibility, providing an efficient and feasible metal artifact removal method for such metal workpieces containing high atomic numbers.

[0006] In order to achieve the above-mentioned purpose, the technical solution adopted by the present invention is as follows: a flowchart of a diffuse metal artifact removal method based on energy spectrum estimation technology is shown in the attached Figure 2 :

[0007] S1: Select a suitable X-ray source according to the specific situation of the test piece, determine the material of the anode target of the X-ray source, and refer to the manual to determine the characteristic peak energy value of the anode material. At the same time, it is necessary to determine the approximate material ratio filled inside the test piece, the structural complexity, size and shape as the prior condition information for the subsequent steps;

[0008] S2: Place the disc-shaped step phantom between the detector and the X-ray source, fix it with a plastic fixing support, and ensure that it is within the imaging field of view. It is necessary to select a suitable material type for the disc-shaped step phantom according to the actual X-ray source energy used. When the X-ray source voltage ranges from 40 KeV to 120 KeV, select aluminum (Al) as the manufacturing material of the disc-shaped step phantom. When the X-ray source voltage ranges from 120 KeV to 300 KeV, select copper (Cu) as the manufacturing material of the disc-shaped step phantom. When the voltage is above 300 KeV, select tungsten (Wu) as the manufacturing material of the disc-shaped step phantom.

[0009] S3: Turn on the X-ray source, obtain the projection data of the step phantom, use the adaptive threshold segmentation method to segment the different grayscales in the projection image, and obtain the "accumulated projection value" of the phantom at different energies to provide a basis for subsequent energy spectrum inference. Including: Assume that the grayscale value range of the image is from I min To I max , then this range needs to be divided into X intervals (where X is generally taken as the energy range of the X-ray source), and the boundary value of each interval can be calculated by the following formula:

[0010]

[0011] Among them, T i is the i∈[1,X-1]th energy threshold, is the width of each energy interval.

[0012] For the grayscale value I(x,y) of each pixel (x,y), it can be divided into different intervals according to these thresholds. The segmentation result can be achieved by the following formula:

[0013]

[0014] Among them, B(x,y) is the interval number after segmentation, and its value range is from 0 to X-1. Indicates a floor operation.

[0015] S4: Use the conjugate gradient method to perform reverse estimation on the segmented projection data B(x, y) obtained in step S3 to obtain the X-ray source energy spectrum curve. Suppose the length of the phantom that the ray that obtains the projection B(x, y) passes through is d k , k is the number of steps of the phantom, the maximum value is equal to the maximum value of the tube voltage of the X-ray source, μ i represents the mass attenuation coefficient of the phantom material at energy i, then:

[0016]

[0017] According to the above formula, we can get i equations containing i unknowns, that is, the actual proportion E of each energy segment of the X-ray source i , the energy spectrum E can be solved by the conjugate gradient method. The above equation can be rewritten as:

[0018] AE=B,

[0019] Among them, A is included The matrix, the elements of which are represented by μ i and d k It is determined that B is a vector containing the measured projection value. Assuming the energy spectrum value of the initial iterative solution is E0, the initial residual sum and search reverse are: r0 = B-AE0, which can be iteratively solved in the following way:

[0020]

[0021] in is the transpose of the residual, p k is the search direction, A is the matrix, and the solution is updated as follows: E k+1 =E k +α k pk ， Update residual method: r k+1 =r k -α.

[0022] S5: remove the phantom and place the workpiece to be tested, scan the workpiece to be tested, and obtain a projection image of the workpiece to be tested;

[0023] S6: According to the energy spectrum data of the X-ray source obtained in step S4, the projection image of the object under test obtained in step S5 is corrected. The projection image needs to be corrected according to the energy spectrum curve of the X-ray source and the prior knowledge of the material filling ratio of the object under test, and finally the beam hardening and photon starvation in the projection image are eliminated to the greatest extent, including the following steps:

[0024] S6.1: Remove Poisson noise and ring artifacts in the projection image of the test piece to prevent errors caused by CT system noise or inconsistent detector response;

[0025] S6.2: Determine the optimal material ratio projection point according to the average gray value of the projection image;

[0026] S6.3: Construct a set of equations based on the prior material ratio and discrete energy spectrum obtained in step S1, solve and use the correction function to correct the projection image, the specific method is as follows:

[0027] Assume that there are Y types of materials in the predicted workpiece, and their proportional function is p i (x), satisfying:

[0028]

[0029] Assuming that all rays are independent and identically distributed, that is, the process of each ray passing through the phantom is independent, and the projection image of each ray has the same distribution, then according to the law of large numbers in statistics, when the number of rays passing through the material is large enough, the pixel average value of the projection image of the phantom can be infinitely close to the projection image of a single ray passing through a workpiece with a real material ratio, that is:

[0030]

[0031] Among them, I avg It is the average pixel value of the projection image of countless rays passing through the phantom, and is calculated as follows:

[0032]

[0033] Then, using this mean as the standard, we search for the closest point in the original projection image as the reference point. We can assume that the material ratio passed by the beam at this point is closest to the actual material ratio of the entire workpiece. Let the material filling ratio of the workpiece being measured at this time be P j , then according to Lambert-Mohr's law, the projection value of the discretized ray passing through the object is:

[0034]

[0035] In the formula, i represents the i-th substance, j represents the j-th energy, X represents the total number of energy spectrum values, Y represents the number of substances contained in the workpiece, and μ ij represents the mass attenuation coefficient of the jth substance at the i-th energy. Assume that the system error shadow can be expressed as:

[0036]

[0037] B(x0,y0) represents the projection value closest to the average value of all projections, and the remaining penetrations are equivalent to the workpiece material filling ratio, which can be solved by the conjugate gradient method again to obtain:

[0038]

[0039] In the formula, B(x i ,y i ) f represents the projection value before correction. Only j+1 equations are needed to solve the proportion p of various substances. j (x), finally, we only need to calculate p j (x) Forward reasoning can get the corrected projection image:

[0040]

[0041] B(x i ,y i ) T Indicates the corrected projected image.

[0042] S7: Reconstruct the corrected projection image of the test piece obtained in step S6 using the FDK algorithm to obtain a high-quality CT image with metal artifacts removed.

[0043] Compared with the traditional metal artifact removal method, the present invention has the following advantages:

[0044] 1) It can solve the diffuse metal artifact problem for which there is no treatment method at present. The method of the present invention removes metal artifacts from the physical mechanism through energy spectrum analysis technology. Compared with the single image post-processing method, it can obtain more information, the processing result is closer to the actual situation, and the image quality is high.

[0045] 2) Strong adaptability. The technology based on energy spectrum estimation can handle metal artifacts of different types and thicknesses and has strong adaptability. Whether it is a large metal object or a small metal fragment, this technology can effectively remove artifacts.

[0046] 3) Simple and efficient. Compared with simple iterative or deep learning metal artifact removal methods, this method uses a phantom to perform reverse estimation of the energy spectrum, without the need for redundant label data or complex iterative reconstruction processes, saving a lot of manual labeling or calculation time, and providing the possibility for CT nondestructive testing of complex workpieces with high atomic numbers. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 It is a flow chart of a diffuse metal artifact removal method based on energy spectrum estimation technology;

[0048] Figure 2 It is a schematic diagram of the disc-shaped stepped phantom structure;

[0049] Figure 3 Schematic diagram of the placement of the disc-shaped stepped phantom in the energy spectrum estimation step. DETAILED DESCRIPTION

[0050] In order to specifically describe the present invention, the technical solution of the present invention is described in detail below in conjunction with the accompanying drawings and specific implementation methods. At the same time, for the purpose of clear and explicit description, implementation examples and corresponding explanations will be given.

[0051] We designed a diffuse metal artifact removal method based on energy spectrum estimation technology. The specific steps are as follows:

[0052] S1: Select a suitable X-ray source according to the specific situation of the test piece, determine the material of the anode target of the X-ray source, and refer to the manual to determine the characteristic peak energy value of the anode material. At the same time, it is necessary to determine the approximate material ratio filled inside the test piece, the structural complexity, size and shape as the prior condition information for the subsequent steps;

[0053] S2: Place the disc-shaped step phantom between the detector and the X-ray source, fix it with a plastic fixing support, and ensure that it is within the imaging field of view. It is necessary to select a suitable material type for the disc-shaped step phantom according to the actual X-ray source energy used. When the X-ray source voltage ranges from 40 KeV to 120 KeV, select aluminum (Al) as the manufacturing material of the disc-shaped step phantom. When the X-ray source voltage ranges from 120 KeV to 300 KeV, select copper (Cu) as the manufacturing material of the disc-shaped step phantom. When the voltage is above 300 KeV, select tungsten (Wu) as the manufacturing material of the disc-shaped step phantom.

[0054] S3: Turn on the X-ray source, obtain the projection data of the step phantom, use the adaptive threshold segmentation method to segment the different grayscales in the projection image, and obtain the "accumulated projection value" of the phantom at different energies to provide a basis for subsequent energy spectrum inference. Including: Assume that the grayscale value range of the image is from I min To I max , then this range needs to be divided into X intervals (where X is generally taken as the energy range of the X-ray source), and the boundary value of each interval can be calculated by the following formula:

[0055]

[0056] Among them, T i is the i∈[1,X-1]th energy threshold, is the width of each energy interval.

[0057] For the grayscale value I(x,y) of each pixel (x,y), it can be divided into different intervals according to these thresholds. The segmentation result can be achieved by the following formula:

[0058]

[0059] Among them, B(x,y) is the interval number after segmentation, and its value range is from 0 to X-1. Indicates a floor operation.

[0060] S4: Use the conjugate gradient method to perform reverse estimation on the segmented projection data B(x, y) obtained in step S3 to obtain the X-ray source energy spectrum curve. Suppose the length of the phantom that the ray that obtains the projection B(x, y) passes through is d k, k is the number of steps of the phantom, the maximum value is equal to the maximum value of the tube voltage of the X-ray source, μ i represents the mass attenuation coefficient of the phantom material at energy i, then:

[0061]

[0062] According to the above formula, we can get i equations containing i unknowns, that is, the actual proportion E of each energy segment of the X-ray source i , the energy spectrum E can be solved by the conjugate gradient method. The above equation can be rewritten as:

[0063] AE=B,

[0064] Among them, A is included The matrix, the elements of which are represented by μ i and d k It is determined that B is a vector containing the measured projection value. Assuming the energy spectrum value of the initial iterative solution is E0, the initial residual sum and search reverse are: r0 = B-AE0, which can be iteratively solved in the following way:

[0065]

[0066] in is the transpose of the residual, p k is the search direction, A is the matrix, and the solution is updated as follows: E k+1 =E k +α k p k , the way to update the residual: r k+1 =r k -α.

[0067] S5: remove the phantom and place the workpiece to be tested, scan the workpiece to be tested, and obtain a projection image of the workpiece to be tested;

[0068] S6: According to the energy spectrum data of the X-ray source obtained in step S4, the projection image of the object under test obtained in step S5 is corrected. The projection image needs to be corrected according to the energy spectrum curve of the X-ray source and the prior knowledge of the material filling ratio of the object under test, and finally the beam hardening and photon starvation in the projection image are eliminated to the greatest extent, including the following steps:

[0069] S6.1: Remove Poisson noise and ring artifacts in the projection image of the test piece to prevent errors caused by CT system noise or inconsistent detector response;

[0070] S6.2: Determine the optimal material ratio projection point according to the average gray value of the projection image;

[0071] S6.3: Construct a set of equations based on the prior material ratio and discrete energy spectrum obtained in step S1, solve and use the correction function to correct the projection image, the specific method is as follows:

[0072] Assume that there are Y types of materials in the predicted workpiece, and their proportional function is p i (x), satisfying:

[0073]

[0074] Assuming that all rays are independent and identically distributed, that is, the process of each ray passing through the phantom is independent, and the projection image of each ray has the same distribution, then according to the law of large numbers in statistics, when the number of rays passing through the material is large enough, the pixel average value of the projection image of the phantom can be infinitely close to the projection image of a single ray passing through a workpiece with a real material ratio, that is:

[0075]

[0076] Among them, I avg It is the average pixel value of the projection image of countless rays passing through the phantom, and is calculated as follows:

[0077]

[0078] Then, using this mean as the standard, we search for the closest point in the original projection image as the reference point. We can assume that the material ratio passed by the beam at this point is closest to the actual material ratio of the entire workpiece. Let the material filling ratio of the workpiece being measured at this time be P j , then according to Lambert-Mohr's law, the projection value of the discretized ray passing through the object is:

[0079]

[0080] In the formula, i represents the i-th substance, j represents the j-th energy, X represents the total number of energy spectrum values, Y represents the number of substances contained in the workpiece, and μ ij represents the mass attenuation coefficient of the jth substance at the i-th energy. Assume that the system error shadow can be expressed as:

[0081]

[0082] B(x0,y0) represents the projection value closest to the average value of all projections, and the remaining penetrations are equivalent to the workpiece material filling ratio, which can be solved by the conjugate gradient method again to obtain:

[0083]

[0084] In the formula, B(x i,y i ) f represents the projection value before correction. Only j+1 equations are needed to solve the proportion p of various substances. j (x), finally, we only need to calculate p j (x) Forward reasoning can get the corrected projection image:

[0085]

[0086] B(x i ,y i ) T Indicates the corrected projected image.

[0087] The present invention is not limited to the above-mentioned specific embodiments, which are merely illustrative and not restrictive. Under the guidance of the present invention, ordinary technicians in this field can make many forms without departing from the scope of protection of the present invention and the claims, all of which are protected by the present invention.

Claims

1. A diffuse metal artifact removal method based on energy spectrum estimation technology, characterized in that: The steps include: S1: Select a suitable X-ray source according to the specific conditions of the test piece, determine the material of the anode target of the X-ray source, refer to the manual to determine the characteristic peak energy value of the anode material, and also determine the approximate proportion of the material filled inside the test piece as the prior condition information for the subsequent steps; S2: Place the disc-shaped step phantom between the detector and the X-ray source, fix it with plastic fixing pillars, and ensure that it is within the imaging field of view; S3: Turn on the X-ray source, obtain the projection data of the step phantom, use the adaptive threshold segmentation method to segment the different grayscales in the projection image, and obtain the "accumulated projection value" of the phantom at different energies to provide a basis for subsequent energy spectrum inference. Including: Assume that the grayscale value range of the image is from I min To I max , then this range needs to be divided into X intervals (where X is generally taken as the energy range of the X-ray source), and the boundary value of each interval can be calculated by the following formula: Among them, T i is the i∈[1,X-1]th energy threshold, is the width of each energy interval. For the grayscale value I(x,y) of each pixel (x,y), it can be divided into different intervals according to these thresholds. The segmentation result can be achieved by the following formula: Among them, B(x,y) is the interval number after segmentation, and its value range is from 0 to X-1. Indicates a floor operation. S4: Use the conjugate gradient method to perform reverse estimation on the segmented projection data B(x, y) obtained in step S3 to obtain the X-ray source energy spectrum curve. Suppose the length of the phantom that the ray that obtains the projection B(x, y) passes through is d k , k is the number of steps of the phantom, the maximum value is equal to the maximum value of the tube voltage of the X-ray source, μ i represents the mass attenuation coefficient of the phantom material at energy i, then: According to the above formula, we can get i equations containing i unknowns, that is, the actual proportion E of each energy segment of the X-ray source i , the energy spectrum E can be solved by the conjugate gradient method. The above equation can be rewritten as: AE=B, Among them, A is included The matrix, the elements of which are represented by μ i and d k It is determined that B is a vector containing the measured projection value. Assuming the energy spectrum value of the initial iterative solution is E0, the initial residual sum and search reverse are: r0 = B-AE0, which can be iteratively solved in the following way: in is the transpose of the residual, p k is the search direction, A is the matrix, and the solution is updated as follows: E k+1 =E k +α k p k , the way to update the residual: r k+1 =r k -α. S5: remove the phantom and place the workpiece to be tested, scan the workpiece to be tested, and obtain a projection image of the workpiece to be tested; S6: Correct the projection image of the test piece obtained in step S5 according to the energy spectrum data of the X-ray source obtained in step S4. S7: Reconstruct the corrected projection image of the test piece obtained in step S6 using the FDK algorithm to obtain a high-quality CT image with metal artifacts removed.

2. The method for removing diffuse metal artifacts based on energy spectrum estimation technology according to claim 1, characterized in that: In step S2, it is necessary to select a suitable material type for the disc-shaped step phantom according to the actual X-ray source energy used. When the X-ray source voltage ranges from 40 Kev to 120 Kev, aluminum (Al) is selected as the manufacturing material of the disc-shaped step phantom. When the X-ray source voltage ranges from 120 Kev to 300 Kev, copper (Cu) is selected as the manufacturing material of the disc-shaped step phantom. When it is above 300 Kev, tungsten (Wu) is selected as the manufacturing material of the disc-shaped step phantom.

3. The method for removing diffuse metal artifacts based on energy spectrum estimation technology according to claim 1, characterized in that: In step S6, the projection image needs to be corrected according to the prior knowledge of the X-ray source energy spectrum curve and the material filling ratio of the object under test, so that the beam hardening and photon starvation phenomena in the final projection image are eliminated to the greatest extent, including the following steps: S6.1: Remove Poisson noise and ring artifacts in the projection image of the test piece to prevent errors caused by CT system noise or inconsistent detector response; S6.2: Determine the optimal material ratio projection point according to the average gray value of the projection image; S6.3: Construct a set of equations based on the prior material ratio and discrete energy spectrum obtained in step S1, solve and use the correction function to correct the projection image, the specific method is as follows: Assume that there are Y types of materials in the predicted workpiece, and their proportional function is p i (x), satisfying: Assuming that all rays are independent and identically distributed, that is, the process of each ray passing through the phantom is independent, and the projection image of each ray has the same distribution, then according to the law of large numbers in statistics, when the number of rays passing through the material is large enough, the pixel average value of the projection image of the phantom can be infinitely close to the projection image of a single ray passing through a workpiece with a real material ratio, that is: Among them, I avg It is the average pixel value of the projection image of countless rays passing through the phantom, and is calculated as follows: Then, using this mean as the standard, we search for the closest point in the original projection image as the reference point. We can assume that the material ratio passed by the beam at this point is closest to the actual material ratio of the entire workpiece. Let the material filling ratio of the workpiece being measured at this time be P j , then according to Lambert-Mohr's law, the projection value of the discretized ray passing through the object is: In the formula, i represents the i-th substance, j represents the j-th energy, X represents the total number of energy spectrum values, Y represents the number of substances contained in the workpiece, and μ ij represents the mass attenuation coefficient of the jth substance at the i-th energy. Assume that the system error shadow can be expressed as: B(x0,y0) represents the projection value closest to the average value of all projections, and the remaining penetrations are equivalent to the workpiece material filling ratio, which can be solved by the conjugate gradient method again to obtain: In the formula, B(x i ,y i ) f represents the projection value before correction. Only j+1 equations are needed to solve the proportion p of various substances. j (x), finally, we only need to calculate p j (x) Forward reasoning can get the corrected projection image: B(x i ,y i ) T Indicates the corrected projected image.

Images