Industrial ct multi-artifact collaborative correction method based on geometric and physical model

Abstract

This invention discloses a collaborative correction method for multiple artifacts in industrial CT based on a geometric and physical model. First, an asymmetric anisotropic spatially variable scattering kernel is constructed based on the geometric parameters and projection gradient information of the industrial CT system. Then, thickness weight data is constructed based on the projection intensity, and an effective attenuation factor is introduced to correct beam hardening. Next, an improved R-L deconvolution iteration is used to couple scattering correction and beam hardening correction. Further, a weighted regularized least squares RWLS model is constructed in the initial reconstruction domain, assigning lower weights to high-density metal regions and their neighborhoods, and combining this with TV penalty solution to correct metal artifacts. Finally, a conventional reconstruction process is performed to obtain the final correction result. This invention incorporates beam hardening artifacts, scattering artifacts, and metal artifacts into a unified geometric and physical model for collaborative processing, solving the problem that individual artifact removal algorithms struggle to eliminate composite artifacts, and improving the grayscale accuracy, edge fidelity, and completeness of composite artifact region correction in industrial CT images.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, specifically relating to a collaborative correction method for multiple artifacts in industrial CT based on geometric and physical models. Background Technology

[0002] Industrial CT often encounters problems such as scattering artifacts, beam hardening artifacts, and metal artifacts when inspecting high-density materials, thick-walled workpieces, complex structural components, and samples containing metal inserts. Among them, scattering artifacts usually appear as fog, decreased contrast, and blurred edges; beam hardening artifacts usually appear as cup-shaped artifacts, dark centers, uneven brightness at the edges, or inconsistent response in thick and thin areas; metal artifacts usually appear as black bands, white bands, radial stripes, and local structural distortion.

[0003] In current industrial CT artifact correction methods, scattering artifacts are typically corrected using empirical symmetric kernels, uniform blur kernels, or simple background estimation methods. Beam hardening artifacts are often corrected using pre-calibrated curves, piecewise fitting, or individual energy spectrum compensation methods. Metal artifacts are mostly corrected using projection interpolation, metal trajectory repair, or reconstruction domain regularization methods. While these methods can improve individual artifacts under specific conditions, in actual industrial CT imaging, scattering artifacts, beam hardening artifacts, and metal artifacts are often not isolated but rather coupled and superimposed, jointly affecting the projection data and reconstruction results. Using only a single artifact correction algorithm is insufficient to completely remove the residual effects in composite artifact regions, especially at high-density edges, areas of abrupt thickness changes, and metal neighborhoods. Often, even after one type of artifact is suppressed, another type may remain or even be amplified.

[0004] Therefore, it is necessary to propose an industrial CT multi-artifact correction method that can jointly handle scattering artifacts, beam hardening artifacts, and metal artifacts. Summary of the Invention

[0005] To address the aforementioned problems, the present invention aims to provide a collaborative correction method for multiple artifacts in industrial CT based on geometric and physical models. This method solves the problems of incomplete correction and significant residual artifacts in composite artifact regions when using only a single artifact algorithm in existing technologies. By constructing an asymmetric anisotropic spatially variable scattering kernel and coupling it with thickness information, effective attenuation information, and uncertainty modeling of metal regions, collaborative correction of scattering artifacts, beam hardening artifacts, and metal artifacts is achieved within a unified geometric and physical model framework. This improves the grayscale accuracy, structural edge fidelity, integrity of composite artifact region correction, and efficiency in engineering applications of industrial CT images.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A collaborative correction method for multiple artifacts in industrial CT based on geometric and physical models includes the following steps: (1) Obtain the original projection data and read the system geometric parameters; (2) Based on the local gradient information of the original projection data and the system geometric parameters, construct an asymmetric anisotropic spatially variable scattering kernel; (3) Construct thickness weight data based on the intensity of the original projection data, and couple the thickness weight data with the asymmetric anisotropic spatial variation scattering kernel to obtain the scattering estimation result after thickness modulation; (4) Introduce an effective attenuation factor to extend the compactness parameter in the scattering estimation result after thickness modulation into a function of the effective attenuation factor, and obtain the beam hardening corrected asymmetric anisotropic spatial variation scattering kernel; (5) Based on the asymmetric anisotropic spatially variable scattering kernel with beam hardening correction, the original projection data is iteratively corrected using an improved RL deconvolution algorithm to obtain scattering and beam hardening coupled corrected projection data. (6) Construct a weighted regularized least squares objective function, and combine it with the total variational penalty function to obtain the final corrected projection data; (7) Reconstruction is performed based on the final corrected projection data to obtain the final corrected industrial CT reconstructed image.

[0007] In this invention, by constructing an asymmetric anisotropic spatially variable scattering kernel and coupling the scattering kernel with thickness information, effective attenuation information, and uncertainties in the metal region, the problem of incomplete correction of composite artifact regions using only a single artifact correction method in the prior art is solved. This invention achieves synergistic correction of scattering artifacts, beam hardening artifacts, and metal artifacts, and the correction results show good grayscale consistency and edge fidelity in composite artifact regions such as high-density edges, abrupt thickness regions, and metal neighborhoods.

[0008] In this invention, the system geometric parameters include the distance from the X-ray source to the detector, the distance from the X-ray source to the rotation center of the sample, the detector pixel size, the projection center coordinates, and the X-ray incident angle.

[0009] In step (2) of this invention, for each pixel (u,v) in the original projection data, a local rotation coordinate system is established with the pixel as the rotation center and the polar angle of the pixel relative to the projection center as the rotation angle. Then, an asymmetric anisotropic space-varying scattering kernel is constructed. ; Where: Norm is the normalization factor, and u and v are the coordinates of the original projected data. , For local rotation coordinates, , For anisotropic scale parameters, α , β For asymmetric direction parameters, orIt is an asymmetric intensity factor. l This is a compactness parameter.

[0010] In this invention, the normalization factor Norm in the asymmetric anisotropic spatially variable scattering kernel is used to ensure that the sum of all elements in the scattering kernel is 1. The specific calculation method is a conventional technique and does not affect the understanding of the technical solution of this invention by those skilled in the art. The asymmetric anisotropic spatially variable scattering kernel model designed in this invention only has a compactness parameter. l While manual settings are required, other asymmetry and orientation information (including anisotropic scale parameters, asymmetric orientation parameters, and asymmetric intensity factors) are automatically determined by the system's geometric parameters and projection data gradients, reducing the difficulty of parameter tuning and the reliance on operator experience.

[0011] In step (3) of this invention, thickness weight data is constructed based on the intensity of the original projection data. By coupling the thickness weight data with an asymmetric anisotropic spatially variable scattering kernel, the scattering estimation results after thickness modulation are obtained. ; Where: I(u,v) is the intensity value of the original projection data at (u,v), I max P represents the maximum intensity value of the original projection data. ie (u,v) represents the estimated projection data of the principal ray at (u,v).

[0012] In this invention, the principal ray projection estimation data P ie (u,v) is obtained from the original projection data after conventional flat-dark field correction and negative logarithmic transformation: , , Where D(u,v) is the dark field projection data at (u,v), that is, the dark current projection data of the detector itself when the X-ray is off and no sample is placed, and F(u,v) is the flat field projection data at (u,v), that is, the air projection data received by the detector when the X-ray is on and no sample is placed. ϵ To prevent small constants from being divided by zero.

[0013] In this invention, thickness weight data W is constructed. s This makes the scattering estimation result positively correlated with the penetration thickness. Specifically, the longer the ray penetration path, the greater the attenuation, and therefore the lower the projection intensity. The corresponding thickness weight value is larger, and the scattering estimation value is larger.

[0014] In step (4), an effective attenuation factor is introduced. The compactness parameter in the thickness-modulated scattering estimation is extended to a function of the effective attenuation factor. A beam-hardened corrected asymmetric anisotropic spatially variable scattering nucleus was obtained. ; Where: I0 is the air intensity value. l base Based on the attenuation factor, c This is the beam hardening correction factor.

[0015] In this invention, the beam hardening correction coefficient is used. c Compactness parameter l Constructed as a function of the effective attenuation factor A, thus setting the compactness parameter. l Convert to set the base attenuation factor l base As the path of the ray penetration increases, A increases. l As the scattering kernel increases, it tends to converge, thus unifying scattering correction and beam hardening correction within the same framework and achieving coupling between the two.

[0016] In step (5), an improved RL deconvolution algorithm is used to iteratively correct the original projection data: a scattering physics model is constructed. According to the formula Iteratively correct the projection data until... Output the current projection data as the projection data for scattering and beam hardening coupling correction; Where: P meas For the original projection data, P ie Principal ray projection estimation data, P n P is the projection data obtained in the nth iteration. n+1 The projection data is obtained in the (n+1)th iteration. ϵ To prevent small constants from being divided by zero, It is an adaptive relaxation factor. , Based on the relaxation factor, Let (u,v) be the gradient magnitude of the original projected data at (u,v). Let P0 be the maximum gradient magnitude of the original projected data, and P0 = P meas。

[0017] In this invention, the improved RL deconvolution algorithm iterative correction adopts a spatial block approximation strategy and utilizes the Fast Fourier Transform (FFT) algorithm combined with parallel computing by a graphics processing unit (GPU) to achieve efficient solution under the condition of spatially variable scattering kernel. Since the scattering kernel S in this invention is spatially variable (i.e., it changes with the spatial position of the projection), it is impossible to directly use the Fast Fourier Transform (FFT) algorithm for convolution acceleration. To address this problem, this invention adopts a spatial block approximation strategy: the projection data is divided into multiple sub-blocks, and within each sub-block, the scattering kernel is approximated as spatially invariant. Based on this, the GPU is used to process each sub-block in parallel, and the Fast Fourier Transform (FFT) algorithm is used within each sub-block to accelerate the convolution calculation. That is, the FFT converts the time-domain convolution into frequency-domain multiplication, significantly reducing the computational complexity and achieving a fast iterative solution for the improved RL deconvolution under the condition of spatially variable scattering kernel.

[0018] For ease of description, the Richardson-Lucy deconvolution algorithm is simply referred to as the RL deconvolution algorithm in this invention, without affecting the understanding of the technical solution of this invention by those skilled in the art. The core of the improved RL deconvolution algorithm in this invention lies in improving the iterative formula of the conventional RL deconvolution algorithm. Specifically, based on the framework of the conventional RL deconvolution algorithm, a new iterative formula is derived according to the scattering physics model of this invention. Furthermore, since conventional RL deconvolution algorithms are prone to overshooting at the edges of high-density samples, specifically manifested as abnormal black stripes (dark areas) appearing immediately adjacent to the edges, this invention introduces an adaptive relaxation factor to prevent overcorrection at the edges of high-density samples. The algorithm adaptively adjusts to the gradient of the projected data: it automatically decreases in high-gradient edge regions to reduce the iteration update step size and suppress black edges caused by overcorrection at high-density edges; it maintains a larger value in flat regions to maintain a normal update step size and quickly remove haze caused by scattering, ultimately yielding the iterative formula. This improvement effectively eliminates black edges and haze while maintaining high-density edge sharpness.

[0019] As is common knowledge, using the RL deconvolution algorithm requires a preset maximum number of iterations, which can be 30, 50, 80, 100, etc. When the number of iterations reaches the maximum number of iterations, the iteration is terminated, and the current projection data is output as the projection data for scattering and beam hardening coupling correction to prevent infinite loops and ensure the real-time performance and reliability of the operation process.

[0020] In step (6), a weighted regularized least squares objective function is constructed. Solve f Obtain the final corrected projection data; in: f For the projection data to be solved, gFor scattering and beam hardening coupling correction of projection data, H is the orthographic projection operator. Here, W is the smoothing parameter, and W is the weighting matrix. Norm g The normalized scattering and beam hardening coupling correction projection data are as follows: TV ( f ) is the total variational penalty function. .

[0021] In this invention, the orthographic projection operator H is calculated using the conventional Siddon algorithm, Joseph algorithm, or variable distance driven algorithm, without affecting the understanding of the technical solution of this invention by those skilled in the art; scattering and beam hardening are coupled together to correct the projection data. g Normalize to the [0,1] interval to obtain normalized scattering and beam hardening coupling correction projection data.

[0022] In this invention, for the sake of formula simplicity, P ideal (u,v) can be written as P ideal S(u) ’ ,v ’ ) can be written as S, S ’ (u ’ ,v ’ ) can be written as S ’ This does not affect the understanding of the technical solution of the present invention by those skilled in the art.

[0023] Furthermore, in step (2), , , ; , , ; , ; ; Where: u c v c Let θ be the coordinates of the projection center, and θ be the rotation angle. c u The horizontal angle of incidence c v The vertical incident angle, n is the incident angle correction factor, SDD is the distance from the X-ray source to the detector, SOD is the distance from the X-ray source to the sample rotation center, Pixel_size is the detector pixel size, and C is a scaling constant related to the energy spectrum. The gradient magnitude is the normalized gradient of the original projected data at (u,v).

[0024] In this invention, the gradient magnitude of the original projection data at (u,v) is... Normalize to the [0,1] interval to obtain the normalized gradient magnitude of the original projected data at (u,v).

[0025] In this invention, the asymmetric direction parameter α , β The rotation angle is adaptively determined based on the local regional characteristics of the projection data: in flat areas, a radial direction (from the projection center to the pixel) is used to simulate the lateral offset of large-angle scattering; in high-density edge areas, a gradient normal direction (perpendicular to the edge) is used to characterize the unilateral scattering tail overflowing from the high-density area to the air side. The rotation angle θ is taken as the polar angle of the current pixel relative to the projection center, so that the anisotropic major axis of the scattering kernel is automatically aligned with the radial geometric direction.

[0026] In this invention, l base = (0.01~10), c =(0.001~1) ϵ =(10 -8 ~10 -4 ), = (0.01~10), = (0.001~1), n= (1~5), C= (0.1~10); More preferably, l base = (0.08~3), c = (0.25~1), ϵ =(10 -8 ~10 -5 ), = (0.08~3), = (0.008~0.3), n= (1~4), C= (0.5~3); More preferably, l base = (0.1~1.5), c =(0.5~1) ϵ =(10 -8 ~10 -6 ), = (0.1~1.5), = (0.01~0.15), n= (1~3), C= (0.8~1.5).

[0027] In step (6) of this invention, the preconditioned conjugate gradient method is used to solve the weighted regularized least squares objective function. As is common knowledge in the field, the preconditioned conjugate gradient method requires the construction of a precondition matrix. In this invention, the precondition matrix of the preconditioned conjugate gradient method is obtained from the image initially reconstructed by the FDK algorithm, and the image initially reconstructed by the FDK algorithm is obtained based on scattering and beam hardening coupling correction projection data. The specific methods for reconstructing the image using the FDK algorithm and for solving the weighted regularized least squares objective function using the preconditioned conjugate gradient method are conventional techniques and do not affect the understanding of the technical solution of this invention by those skilled in the art.

[0028] In this invention, the original projection data is downsampled, and steps (2) to (6) are performed on the downsampled projection data to obtain low-resolution corrected projection data. The low-resolution corrected projection data is upsampled and mapped to the original projection data to obtain the final corrected projection data. Then, step (7) is performed to obtain the final corrected industrial CT reconstructed image.

[0029] This invention discloses a storage medium storing a processor-executable program for executing the above-described industrial CT multi-artifact collaborative correction method based on geometric and physical models.

[0030] The present invention discloses an electronic device, including a memory and a processor, characterized in that: the memory stores a method for executing the above-mentioned industrial CT multi-artifact collaborative correction method based on geometric and physical models.

[0031] This invention discloses the application of the industrial CT image artifact correction method in correcting image artifacts.

[0032] Due to the application of the above technical solution, the beneficial effects of the present invention are as follows: (1) Based on geometric and physical models, this invention creatively proposes a composite artifact correction method. Compared with the existing technology that only uses a single artifact correction algorithm, it can effectively reduce residual interference in the composite artifact region and improve the integrity and consistency of artifact removal. It is especially suitable for projections with composite artifact problems. (2) The scattering kernel designed in this invention is not a simple conventional symmetrical fixed kernel, but a scattering kernel with asymmetry, anisotropy and spatial variation, which is more in line with the real scattering distribution law under the combined effect of oblique incidence, cone-beam geometry, thickness gradient and boundary changes in industrial CT, and can improve the physical rationality and correction accuracy of scattering estimation. (3) This invention combines the scattering kernel with the thickness weight and couples the scattering correction with the beam hardening correction through an effective attenuation factor, so that the scattering kernel not only reflects the spatial diffusion morphology, but also reflects the scattering change characteristics under different penetration thicknesses and energy spectrum hardening states, thereby improving the correction consistency under thick-walled, high-density sample conditions. (4) This invention uses a weighted regularized least squares (RWLS) objective function, a weighted matrix and a total variation (TV) penalty to adaptively weaken unreliable projection data in the metal region and its neighborhood during the reconstruction process. This can effectively suppress shadows, stripes and local distortions near the metal edge and improve the reconstruction quality of the metal neighborhood. (5) Furthermore, the present invention adopts an engineering implementation method of downsampling pre-calculation and correction data upsampling backfilling, which significantly reduces the computational burden of high-resolution direct iteration while ensuring the correction effect, and is suitable for industrial field applications. Attached Figure Description

[0033] Figure 1 This is a flowchart of the method steps of the present invention.

[0034] Figure 2 This is a comparison chart showing the calibration effects of alloy blade samples. Among them: Figure 2 (a) is the original reconstructed image. Figure 2 (b) is the reconstructed image after correction of beam hardening artifacts only. Figure 2 (c) is the reconstructed image after correction of only scattering artifacts. Figure 2 (d) is the reconstructed image after correction of only metal artifacts. Figure 2 (e) is the reconstructed image after correction using the method of the present invention. Detailed Implementation

[0035] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0036] The spatial block approximation strategy, Fast Fourier Transform (FFT) algorithm, graphics processing unit (GPU), calculation method of orthographic projection operator H, FDK algorithm, preconditional conjugate gradient method (PCG), and upsampling and downsampling methods used in this invention are all prior art and do not affect the understanding of this invention by those skilled in the art. The inventiveness of this invention lies in its collaborative correction framework, constructed through unified scattering kernel modeling, thickness and effective attenuation coupling, and weighted regularization constraints. This framework solves the problem of incomplete sequential stitching correction of composite artifact regions in existing technologies, reduces the difficulty of parameter tuning and the reliance on operator experience, and is suitable for industrial field applications. Example 1

[0037] A collaborative correction method for multiple artifacts in industrial CT based on geometric and physical models includes the following steps: (1) Obtain the original projection data and read the system geometric parameters. Before acquiring the raw projection data, the industrial CT system is first calibrated and its geometric parameters are determined. The specific procedures are standard techniques, as shown below: ①System calibration: Dark current correction involves acquiring dark current images without radiation exposure and subtracting them from the original projection. Gain correction involves acquiring a flat-field image of uniformly irradiated light and normalizing the original projection to eliminate detector response inhomogeneity. Bad pixel and bad column correction: Repairing failed detector pixels through neighborhood interpolation; ② Geometric parameter calibration: distance from X-ray source to detector (SDD), distance from X-ray source to sample rotation center (SOD), detector pixel size (Pixel_size), detector dimensions (a×b, where a is the number of detector pixel rows and b is the number of detector pixel columns), rotation axis tilt correction, and detector in-plane tilt correction; ③ Calibration of the energy spectrum characteristics of the X-ray source, including calibration of the effective energy spectrum distribution and focal spot size; After system calibration, the calibrated original projection data is acquired. The original projection data is then routinely downsampled, and system geometric parameters are read, including the distance from the X-ray source to the detector (SDD), the distance from the X-ray source to the sample rotation center (SOD), the detector pixel size (Pixel_size), and the projection center coordinates (u...). c ,v c ), u c =a / 2、v c =b / 2), angle of incidence of ray ( , This is used to determine the size, orientation, and spatial distribution of the scattering nuclei in subsequent steps.

[0038] (2) Constructing asymmetric anisotropic space-varying scattering kernels ① Calculate the basic scale parameters: ; Where: C is a proportionality constant, C=0.8; ②Based on the horizontal angle of incidence c u Vertical angle of incidence c v Calculate the anisotropic scale parameters: , ; Where: n is the incident angle correction factor, n=2; ③ Perform gradient calculation on the original projection data to obtain the gradient magnitude. The gradient magnitude of the original projected data at (u,v) is obtained by normalizing it to the [0,1] interval, thus yielding the normalized gradient magnitude of the original projected data. Used to distinguish between flat areas and edge areas: in flat areas ( or Approaching 0), using the radial direction; in the edge region ( or(If the direction is relatively large), switch to the gradient normal direction to obtain the asymmetric direction parameters: , ; Where u and v are the coordinates of the original projection data. Let be the normalized gradient magnitude of the original projected data at (u,v). ; ④ For each pixel (u,v) in the original projection data, rotate around that pixel as the center of rotation and its relative projection center (u c ,v c Using the polar angle θ as the rotation angle, establish a local rotating coordinate system. Constructing asymmetric anisotropic space-varying scattering kernels: ; in: , , , Norm is a normalization factor used to ensure that the sum of all elements in the scattering kernel K is 1. l This is a compactness parameter.

[0039] (3) Thickness direction modulation ① Calculate the thickness weight data based on the intensity I of the original projection data: ; Where: I(u,v) is the intensity value of the original projection data at (u,v), I max The maximum intensity value of the original projection data; ② Multiply the thickness-weighted data with the ideal projection data, and then convolve it with the asymmetric anisotropic spatially variable scattering kernel to obtain the thickness-modulated scattering estimation result: ; Where: P ie (u,v) represents the estimated projection data of the principal ray at (u,v), obtained from the original projection data after conventional flat-field correction and negative logarithmic transformation. The specific steps are as follows: First, perform flat-field correction to obtain the projection data after flat-field correction. ; Then, a negative logarithmic transformation is performed to obtain the estimated data for the principal ray projection. ; D(u,v) represents the dark field projection data at (u,v), which is the dark current projection data of the detector itself when the X-ray is off and no sample is placed. F(u,v) represents the flat field projection data at (u,v), which is the air projection data received by the detector when the X-ray is on and no sample is placed. ϵ To prevent small constants from being divided by zero, ϵ =10 -6 .

[0040] (4) Beam hardening correction Introducing an effective attenuation factor The compactness parameter in the thickness-modulated scattering estimation result is set as a function of the effective attenuation factor: ; Thus, a beam-hardened corrected asymmetric anisotropic space-varying scattering nucleus is obtained: ; Where: I0 is the air intensity value. l base Based on the attenuation factor, l base =0.15, c This is the beam hardening correction factor. c =1. The greater the attenuation of the ray as it penetrates deeper into the path, the larger A becomes. l The larger the value, the more convergent the scattering nucleus (i.e., the more concentrated the scattered energy is in the forward direction).

[0041] (5) Scattering deconvolution correction Based on a beam-hardening corrected asymmetric anisotropic spatially varying scattering kernel, an improved RL deconvolution algorithm is used to process the original projection data P. meas Perform iterative correction: Constructing a scattering physics model: Using the original projection data as the initial value for iteration, P0=P meas The projected data is updated in each iteration according to the following formula: until The output satisfies P at time n+1 As scattering and beam hardening coupling correction projection data g; Where: P meas For the original projection data, P ie Principal ray projection estimation data, P n P is the projection data obtained in the nth iteration. n+1 The projection data is obtained in the (n+1)th iteration. ϵ To prevent small constants from being divided by zero, ϵ =10 -6 , It is an adaptive relaxation factor. , Based on the relaxation factor, =0.15, Let (u,v) be the gradient magnitude of the original projected data at (u,v). The maximum gradient magnitude of the original projection data; In this invention, the improved RL deconvolution algorithm iterative correction adopts a spatial block approximation strategy, which divides the projected data into multiple local homogeneous regions. In each region, it approximates a spatially invariant kernel, uses the Fast Fourier Transform (FFT) algorithm to accelerate the convolution calculation, and combines a graphics processing unit (GPU) to perform parallel calculations on each region to improve the solution efficiency.

[0042] (6) Weighted correction of high-density areas of metal artifacts ① The FDK algorithm was used to initially reconstruct the image from the scattering and beam hardening coupling correction projection data; ② Construct the weighted regularized least squares (RWLS) objective function: ; in: f Let g be the projection data to be solved, and g be the projection data corrected by scattering and beam hardening coupling. For smoothing parameters, =0.01, W is the weighting matrix, , TV ( f ) is the total variational penalty function. , H is the orthographic projection operator, calculated using the conventional Siddon algorithm, Joseph algorithm, or variable distance-driven algorithm, without affecting the understanding of the technical solution of this invention by those skilled in the art. Correcting projection data by coupling scattering with beam hardening g Normalize to the [0,1] interval to obtain normalized scattering and beam hardening coupling correction projection data; ③ Solve the weighted regularized least squares objective function using the preconditioned conjugate gradient method (PCG). J ( f ): Construct a preconditioning matrix M to approximate the feature value distribution of the image initially reconstructed by the FDK algorithm, and then perform variable substitution. f = Mf' The original objective function is rewritten as follows: Simultaneously, smoothing parameters are used in high-density metal regions. Reduce TV penalty intensity to avoid excessive smoothing of metal edges due to global TV penalty; This strategy can improve the matrix condition number and accelerate the iterative solution of PCG. The solution yields... f’ Then, according to f = Mf' Obtain projection data f .

[0043] (7) Reconstruction The projection data is conventionally upsampled to the original high resolution and mapped back to the original projection data to obtain the final corrected projection data. Conventional reconstruction is then performed based on the final corrected projection data to obtain the final corrected industrial CT reconstructed image. As is common practice, using the RL deconvolution algorithm requires a preset maximum number of iterations. In this embodiment, the maximum number of iterations for the RL deconvolution algorithm is 30. When the number of iterations reaches the maximum number of iterations, the iteration is terminated, and the current projection data is output as the projection data for scattering and beam hardening coupling correction to prevent infinite loops and ensure the real-time performance and reliability of the operation process. Example 2

[0044] A storage medium storing a processor-executable program for executing the industrial CT multi-artifact collaborative correction method based on a geometric and physical model according to Embodiment 1. Example 3

[0045] An electronic device includes a memory and a processor, the memory storing a method for performing collaborative correction of multiple artifacts in industrial CT based on a geometric and physical model, as described in Embodiment 1.

[0046] Comparative Example 1 The original projection was corrected using a conventional energy spectrum compensation method to correct beam hardening artifacts.

[0047] Comparative Example 2 The original projection is corrected for scattering artifacts using a conventional empirical symmetric kernel method.

[0048] Comparative Example 3 The original projection was repaired and corrected using a conventional metal trajectory method to correct metal artifacts.

[0049] Application Examples Alloy blades (which have cavity structures, local thickness variation areas, and edge transition areas, reflecting the superposition effects of scattering artifacts, beam hardening artifacts, and local stripe artifacts in actual samples) were used as test samples for conventional industrial CT scanning and reconstruction: Scanning equipment: XTremeVista 3000 device from Micro-Kuang Technology (Suzhou) Co., Ltd.; Scanning parameters: tube voltage 440kV, tube current 110μA, exposure time 1000ms, number of projection frames 1600, filter is 10mm Cu, reconstructed voxel size is 55.469μm; Geometric parameters: the distance from the X-ray source to the detector (SDD) is 857.14 mm, the distance from the X-ray source to the center of rotation of the sample (SOD) is 342.05 mm, the detector pixel size (Pixel_size) is 139 μm, the detector size is 3032×3032, and the projection center coordinates are (1516, 1516). like Figure 2 (a) shows the original reconstructed image of the alloy blade sample. It can be seen that there are obvious gray fog, uneven gray levels, and local artifact residues in the thick-walled region, edge transition region, and internal cavity neighborhood of the alloy blade.

[0050] The original projection data were processed using the correction methods of Example 1, Comparative Example 1, Comparative Example 2, and Comparative Example 3, respectively. The image reconstructed using Example 1 is shown below. Figure 2 (e) See the comparative one-correction reconstructed image. Figure 2 (b) See the comparative two-correction reconstructed image. Figure 2 (c) See the comparative three-correction reconstructed image. Figure 2 (d). It can be seen that after correction using the method of the present invention, the fog phenomenon in the relevant area is significantly reduced, the local gray-scale consistency is improved, the edge transition is clearer, and the stripes and local artifacts are suppressed. However, after using the separate correction method, although some artifacts are reduced to a certain extent, there are still residual artifacts, gray-scale inconsistencies, or insufficient local correction in local high-density edge areas and areas with abrupt changes in thickness.

[0051] In theory, a combination of conventional energy spectrum compensation correction, empirical symmetry kernel correction, and metal trajectory repair correction can be performed on the original projection. However, in practice, it has been found that the correction results of the previous stage will affect the correction effect of the next stage, resulting in mutual constraints between the corrections. This leads to high engineering implementation complexity and makes it difficult to deploy and stably apply the corrections in industrial sites.

[0052] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for industrial CT multi-artifact collaborative correction based on geometric and physical models, characterized in that, Includes the following steps: (1) Obtain the original projection data and read the system geometric parameters; (2) Based on the local gradient information of the original projection data and the system geometric parameters, construct an asymmetric anisotropic spatially variable scattering kernel; (3) Construct thickness weight data based on the intensity of the original projection data, and couple the thickness weight data with the asymmetric anisotropic spatial variation scattering kernel to obtain the scattering estimation result after thickness modulation; (4) Introduce an effective attenuation factor to extend the compactness parameter in the scattering estimation result after thickness modulation into a function of the effective attenuation factor, and obtain the beam hardening corrected asymmetric anisotropic spatial variation scattering kernel; (5) Based on the asymmetric anisotropic spatially variable scattering kernel with beam hardening correction, the original projection data is iteratively corrected using an improved RL deconvolution algorithm to obtain scattering and beam hardening coupled corrected projection data. (6) Construct a weighted regularized least squares objective function, and combine it with the total variational penalty function to obtain the final corrected projection data; (7) Reconstruction is performed based on the final corrected projection data to obtain the final corrected industrial CT reconstructed image.

2. The method of claim 1, wherein the method further comprises: The system's geometric parameters include the distance from the X-ray source to the detector, the distance from the X-ray source to the sample's rotation center, the detector's pixel size, the projection center coordinates, and the X-ray incident angle.

3. The industrial CT multi-artifact collaborative correction method based on geometric and physical models according to claim 1, characterized in that: In step (2), for each pixel (u,v) in the original projection data, a local rotation coordinate system is established with the pixel as the rotation center and the polar angle of the pixel relative to the projection center as the rotation angle. Then, an asymmetric anisotropic space-varying scattering kernel is constructed. ; Where: Norm is the normalization factor, and u and v are the coordinates of the original projected data. , For local rotation coordinates, , For anisotropic scale parameters, α , β For asymmetric direction parameters, η It is an asymmetric intensity factor. λ For compactness parameters; In step (3), thickness weight data is constructed based on the intensity of the original projection data. By coupling the thickness weight data with an asymmetric anisotropic spatially variable scattering kernel, the scattering estimation results after thickness modulation are obtained. ; where: I(u, v) is the intensity value of the original projection data at (u, v), I max is the maximum intensity value of the original projection data, P ie (u, v) is the primary ray projection estimation data at (u, v). In step (4), an effective attenuation factor is introduced The compactness parameter in the thickness-modulated scatter estimate is extended to a function of the effective attenuation factor , resulting in a beam-hardening-corrected asymmetric anisotropic scatter kernel ; Where: I0 is the air intensity value. λ base Based on the attenuation factor, γ This is the beam hardening correction factor; In step (5), an improved RL deconvolution algorithm is used to iteratively correct the original projection data: a scattering physics model is constructed. According to the formula Iteratively correct the projection data until... Output the current projection data as the projection data for scattering and beam hardening coupling correction; Where: P meas For the original projection data, P ie Principal ray projection estimation data, P n P is the projection data obtained in the nth iteration. n+1 The projection data is obtained in the (n+1)th iteration. ϵ To prevent small constants from being divided by zero, As an adaptive relaxation factor, , Based on the relaxation factor, Let (u,v) be the gradient magnitude of the original projected data at (u,v). Let P0 be the maximum gradient magnitude of the original projected data, and P0 = P meas ; In step (6), a weighted regularized least squares objective function is constructed. Solve f Obtain the final corrected projection data; in: f For the projection data to be solved, g For scattering and beam hardening coupling correction of projection data, H is the orthographic projection operator. Here, W is the smoothing parameter, and W is the weighting matrix. Norm g The normalized scattering and beam hardening coupling correction projection data are as follows: TV ( f ) is the total variational penalty function. .

4. The method of claim 3, wherein the method further comprises: In step (2), ， ， ； ， ， ； ， ； ； Where: u c v c Let θ be the coordinates of the projection center, and θ be the rotation angle. γ u The horizontal angle of incidence γ v The vertical incident angle, n is the incident angle correction factor, SDD is the distance from the X-ray source to the detector, SOD is the distance from the X-ray source to the sample rotation center, Pixel_size is the detector pixel size, and C is a scaling constant related to the energy spectrum. The gradient magnitude is the normalized gradient of the original projected data at (u,v).

5. The industrial CT multi-artifact collaborative correction method based on geometric and physical models according to claim 1, characterized in that: In step (6), the preconditioned conjugate gradient method is used to solve the weighted regularized least squares objective function.

6. The method of claim 1, wherein: The precondition matrix of the precondition conjugate gradient method is obtained from the image initially reconstructed by the FDK algorithm.

7. The method of claim 1, wherein: The original projection data is downsampled, and steps (2) to (6) are performed on the downsampled projection data to obtain low-resolution corrected projection data. The low-resolution corrected projection data is upsampled and mapped to the original projection data to obtain the final corrected projection data. Then, step (7) is performed to obtain the final corrected industrial CT reconstructed image.

8. A storage medium storing a program executable by a processor, characterized by: The program is used to execute the industrial CT multi-artifact collaborative correction method based on geometric and physical models as described in any one of claims 1-7.

9. An electronic device comprising a memory, a processor, characterized in that: The memory stores the method for performing the industrial CT multi-artifact collaborative correction method based on a geometric and physical model as described in any one of claims 1-7.

10. The application of the industrial CT image artifact correction method of claim 1 in correcting image artifacts.

Images