Unsupervised inter-layer complementarity enhanced ring artifact removal method

Abstract

The invention discloses an unsupervised method for removing ring artifacts with enhanced interlayer complementarity. According to the method, a polynomial coefficient compensation equation is defined to optimize the non-uniform response of the detector based on the principle that ring artifacts are generated based on the non-uniform response of the detector. In consideration of randomness of detector response and independence in the vertical direction, a 2.5 D data input format is designed, remarkable interlayer information between adjacent slices is utilized, interlayer information constraint is realized by applying a TV minimization algorithm, and preliminary ring removal and denoising are carried out. In the polar coordinate system, the ring artifact is represented as a vertical stripe, the gradient of the stripe in the horizontal direction is large, a horizontal differential operator is designed to extract the polar coordinate image gradient, the loss is calculated, and the compensation coefficient is updated. According to the method, hyper-parameters do not need to be adjusted, images do not need to be trained, the ring and noise removal effect is good, and the automatic process is completely achieved.

Description

Technical Field

[0001] The invention designs an unsupervised inter-layer complementarity enhanced ring artifact removal method, which belongs to the field of computer image processing. Background Art

[0002] Computed tomography systems equipped with photon counting detectors have shown great potential in the field of image reconstruction. However, the inconsistency of detector response leads to stripe artifacts in the projection data, which in turn causes concentric ring artifacts in the image, seriously affecting the image quality and texture details.

[0003] There are two types of ring artifacts: high-frequency ring artifacts and low-frequency ring artifacts. High-frequency ring artifacts appear as dense concentric circles, while low-frequency ring artifacts appear as band-shaped circular structures with a certain width. High-frequency ring artifacts are mainly caused by three factors: detector dead spots, non-uniform response between detector elements, and gaps between detectors. Low-frequency ring artifacts are usually caused by discontinuities between different pixel clusters of the detector.

[0004] Traditional ring artifact removal work can be roughly divided into hardware-based and software-based methods. Hardware-based methods change the relative position of the detector and the scanned object to offset the non-uniformity of the detector response. However, these methods require additional mechanisms to achieve relative motion, increase the manufacturing difficulty and cost of the CT system, and may amplify the noise of the reconstructed CT image. Software methods are divided into pre-processing methods, post-processing methods, and optimization-based methods. The pre-processing method is to remove vertical stripe artifacts in the sinusoidal image so that the reconstructed image no longer has ring artifacts. The post-processing method is to transform the image with ring artifacts from the Cartesian domain to the polar coordinate domain so that the artifact is displayed in the form of stripes, which is easy to detect and eliminate. The optimization-based method constructs an objective function according to the nature of the artifact and obtains the corrected image by minimizing the function. Regardless of which method is used, it requires rich knowledge and experience in image reconstruction and filtering, and needs to adjust hyperparameters according to image data, which is not suitable for promotion as a universal method. Therefore, the present invention proposes an automated ring artifact removal method based on detector non-uniform response compensation. The method does not require any pre-data and training hyperparameters, and while achieving ring removal, it also has a certain inhibitory effect on reconstruction noise. Summary of the invention

[0005] In view of the problems of cumbersome steps and difficult hyperparameter adjustment in traditional ring removal methods, this paper proposes a ring removal method based on detector non-uniform response compensation. The ring artifacts appear as vertical stripes in the polar coordinate system, and the horizontal differential Sobel operator and Scharr operator are used to extract the stripe gradient to calculate the loss function to optimize the learnable detector compensation coefficient equation. The algorithm calculates the minimum loss value in an iterative manner and continuously compensates for the non-uniform response of the detector.

[0006] The present invention uses polynomial compensation to realize the non-uniform compensation of the detector pixel response, and takes TV as an example (not limited to TV) to realize the inter-layer information constraint update. Assuming that the maximum number of iterations of the inter-layer information constraint update is K, and the maximum number of iterations of the polynomial compensation update is M, the present invention describes three inter-layer information constraint update and polynomial compensation update combination strategies, such as Figure 1 As shown. They are:

[0007] (1) Set the polynomial compensation update to external iteration and the inter-layer information constraint update to internal iteration. The polynomial compensation is iterated M times, and the inter-layer information constraint is iterated K·M times. The advantage of this combination is that K can be set to a small value, which can suppress the distortion problem of smoothed images while ensuring a good overall denoising effect.

[0008] (2) The idea of ​​​​series is adopted. The polynomial compensation is iterated M times first, and the inter-layer information is updated and iterated K times. However, this combination does not take into account the distortion problem caused by the TV minimization algorithm. If K is too large, it may cause image distortion; if K is too small, the TV minimization algorithm has poor de-ringing and denoising effects.

[0009] (3) Set the number of external iterations T. Both polynomial compensation iterations and inter-layer information constraint updates are internal iterations and are updated alternately. For one external iteration, first perform M polynomial compensation iterations, then perform K inter-layer information constraint updates. A total of T·M polynomial compensation iterations and T·K inter-layer information constraint iterations are performed.

[0010] In order to achieve the above object, the technical solution of the present invention is as follows: an unsupervised ring artifact removal method with enhanced inter-layer complementarity. The method comprises the following steps:

[0011] Step 1: Design a 2.5D input format based on inter-slice complementarity. The detector responses are independent in the head-to-foot direction, i.e., the z-axis direction. The detector responses are random. This non-uniform response means that the annular artifact information of each slice is complementary. This method proposes a 2.5D data input structure, takes a Z-axis slice group (the present invention is a group of 5 slices), and places the number of slices in the channel dimension. Compared with traditional 2D and 3D structures, it improves computational efficiency while retaining effective spatial redundant information. Step 2: Define a polynomial compensation coefficient equation. The main cause of the ring artifact is the non-uniform response between the detector elements and the detector pixel blocks. The compensation coefficient polynomial is defined to optimize the non-uniform response of the detector and realize the correction of the ring artifact. For the j-th pixel of the detector, use formula (1) for optimization.

[0012]

[0013] Among them, S j is the original pixel value of the jth pixel of the detector, S j,corr is the pixel value after compensation, a j,i is the i-th order compensation coefficient of the j-th pixel of the detector. This parameter is a differentiable independent variable parameter with gradient. Initialize a j,1 =1, and 0 for the others. Here i=3. Taking the experimental data in this method as an example, the size of the sinusoidal projection image is 600x384x59. That is, there are 59 slices, 600 projection angles, and 384 detector pixels. Since the response value of the detector is the same at each angle, the parameters that need to be learned to process an image slice are only 384x4.

[0014] Step 3: Use the Total Variation (TV) minimization method to introduce inter-slice complementary information and suppress noise. The introduction of total variation information in the Z-axis direction can combine the redundant information of artifacts in different slices to achieve the effect of removing some ring artifacts. The total variation is the sum of the differences in adjacent pixel values ​​in an image, which reflects the smoothness or edge information of the image. By minimizing the total variation of the image, noise and redundant information can be suppressed while retaining the important features of the image, as shown in Equations 2 and 3.

[0015]

[0016] img corr =img+Norm(J) (3)

[0017] Among them, the first is the total variation of the image. Limiting the total variation can suppress noise. The second term is the fidelity term, which prevents the denoised image from being too different from the original image. λ is the relaxation factor. Since K is small, the image distortion phenomenon is not obvious, so it is set to λ=0 in this algorithm.

[0018] Step 4: Use Gaussian strided convolution to obtain multi-scale feature maps. There are two types of ring artifacts: high-frequency rings and low-frequency rings. High-frequency ring artifacts are densely distributed, while low-frequency rings are more dispersed and have a certain width. In previous methods, low-frequency rings are often difficult to remove. This method proposes Gaussian strided convolution to achieve downsampling. Since the input data dimension is [B, S, H, W], in order to ensure that each slice is independent of each other during downsampling, group convolution is used. Use the Gaussian kernel as the convolution kernel and set the convolution Stride = 2 to achieve image downsampling. Through strided convolution, image smoothing can be achieved, and some low-frequency rings can be converted into high-frequency rings to achieve better de-ringing effect.

[0019] Step 5: Propose an online dual-domain conversion algorithm. The present invention proposes an online dual-domain conversion algorithm that can realize the conversion between the Cartesian domain and the polar coordinate domain while retaining the gradient information.

[0020] Cartesian2Polar: Defines the angle θ and radius r in the polar coordinate system. Use torch.meshgrid to combine θ and r into a network that includes the coordinates of each point in polar coordinates. Use the torch.polar function to convert the grid coordinates to polar coordinates, as shown in Equation 3. Finally, use the F.grid_sample method to sample the Cartesian domain image to the polar coordinate domain.

[0021] z=r·cos(θ)+r·sin(θ)·j (4)

[0022] Where r is the radius and θ is the angle (radians).

[0023] Polar2Cartesian: Since the origin of a polar coordinate image is the center of the image, the image is split into a left image and a right image along the central axis, the left image is flipped horizontally, and the two images are spliced ​​vertically. Create a grid, calculate the Cartesian coordinates (X, Y) corresponding to each pixel, and convert the X and Y coordinates to coordinates relative to the center of the image. Calculate the distance r and angle θ from each pixel to the center of the image. Use the F.grid_sample function to sample the spliced ​​image according to the calculated polar coordinate network, obtain the converted image, and adjust the image.

[0024] Step 6: Use the horizontal differential operator to extract the stripe artifacts in the polar coordinate domain and calculate the loss function. The loss function generation process is shown in Figure 8After the artifact images at different resolution scales are converted to polar coordinate images through the online coordinate domain conversion algorithm, the horizontal gradient is calculated using the horizontal differential operator. Since the gradient at the stripe artifact is large, the loss function can capture the stripe artifact well. By minimizing the loss function, the stripe artifact is removed in the polar coordinate domain, that is, the ring artifact is removed in the image domain.

[0025] (4) Step 7: Iterate the above process 50 times, the ring artifacts and noise are effectively suppressed, and a ring-free high-quality image is output.

[0026] Compared with the prior art, the advantages of the present invention are as follows: (1) Good de-ringing and denoising effects. The present invention proposes an unsupervised ring artifact removal method with enhanced inter-layer complementarity. Starting from the principle of ring artifact generation, the detector response is optimized by polynomial compensation coefficients, and the de-ringing effect is good. At the same time, the TV algorithm is combined in the iterative process to make full use of inter-chip information for de-ringing, which has a certain denoising effect. (2) It is universal and parameter-free. This method does not require the adjustment of any hyperparameters, does not require any professional knowledge and experience, and has a good artifact removal effect for any detector device. (3) It runs fast. The method runs based on GPU, has few learning parameters, occupies little memory, has the potential to run on low-computing power GPU devices, and requires a short time for actual testing. (4) The present invention is an automated, end-to-end method that can be integrated into the reconstruction software of the CT scanner. Only the corresponding button needs to be designed to obtain high-quality images with ring artifacts removed. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 Polynomial compensation update and inter-layer information constraint update compensation combination strategy diagram

[0028] Figure 2 Pseudo code for the overall algorithm;

[0029] Figure 3 is a flow chart of the present invention;

[0030] Figure 4 Update pseudocode for inter-layer information constraints;

[0031] Figure 5 is the complementarity diagram between slices; (a) is the non-uniform response of the detector, (b) to (d) are the extracted continuous slices;

[0032] Figure 6 is the relationship diagram between the detector non-uniform response and the ring artifact;

[0033] Figure 7 This is the effect diagram of the online dual-domain conversion algorithm;

[0034] Figure 8 Schematic diagram of Gaussian convolution kernel pooling;

[0035] Fig. 9 Schematic diagram for loss function calculation;

[0036] Fig.10 The test results of the algorithm. DETAILED DESCRIPTION

[0037] The present invention is further explained below in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention. After reading the present invention, various equivalent forms of modifications to the present invention by those skilled in the art all fall within the scope defined by the claims attached to this application.

[0038] Example 1: An unsupervised ring artifact removal method with enhanced inter-layer complementarity, a combined strategy of updating polynomial coefficients and updating inter-layer information constraints is shown in Figure 1 , the following steps adopt scheme (1). The detailed algorithm is shown in Figures 2 to 4 The method comprises the following steps:

[0039] Step 1: Design a 2.5D input format based on inter-chip complementarity. The detector responses are independent in the head-foot direction, i.e., the z-axis direction. The detector responses are random, such as Figure 5 This non-uniform response means that the annular artifact information of each slice is complementary. This method proposes a 2.5D data input structure, takes a Z-axis slice group (the present invention is a group of 5 slices), and places the number of slices into the channel dimension. Compared with traditional 2D and 3D structures, it improves computational efficiency while retaining effective spatial redundant information.

[0040] Step 2: Define the polynomial compensation coefficient equation. The main reason for the ring artifact is the non-uniform response between detector elements and detector pixel blocks, see Figure 6 The compensation coefficient polynomial is defined to optimize the non-uniform response of the detector and realize the correction of the ring artifact. For the j-th pixel of the detector, the optimization is performed using equation (1).

[0041]

[0042] Among them, S j is the original pixel value of the jth pixel of the detector, S j,corr is the pixel value after compensation, a j,i is the i-th order compensation coefficient of the j-th pixel of the detector. This parameter is a differentiable independent variable parameter with gradient. Initialize a j,1=1, and 0 for the others. Here i=3. Taking the experimental data in this method as an example, the size of the sinusoidal projection image is 600x384x59. That is, there are 59 slices, 600 projection angles, and 384 detector pixels. Since the response value of the detector is the same at each angle, the parameters that need to be learned to process an image slice are only 384x4.

[0043] Step 3: Use the Total Variation (TV) minimization method to introduce inter-slice complementary information and suppress noise. Introduce the total variation information in the Z-axis direction, which can combine the redundant information of artifacts in different slices to achieve the effect of removing some ring artifacts. The total variation is the sum of the differences in adjacent pixel values ​​in the image, which reflects the smoothness or edge information of the image. By minimizing the total variation of the image, noise and redundant information can be suppressed while retaining the important features of the image, as shown in Equations 2 and 3.

[0044]

[0045] img corr =img+Norm(J) (3)

[0046] Among them, the first is the total variation of the image. Limiting the total variation can suppress noise. The second term is the fidelity term, which prevents the denoised image from being too different from the original image. λ is the relaxation factor. Since K is small, the image distortion phenomenon is not obvious, so it is set to λ=0 in this algorithm.

[0047] Step 4: Use Gaussian strided convolution to obtain multi-scale feature maps. Ring artifacts are divided into high-frequency rings and low-frequency rings. High-frequency ring artifacts are densely distributed, while low-frequency rings are more dispersed and have a certain width. In previous methods, low-frequency rings are often difficult to remove. This method proposes Gaussian strided convolution to achieve downsampling, such as Figure 7 . Since the input data dimension is [B, S, H, W], group convolution is used to ensure that each slice is independent during downsampling. The Gaussian kernel is used as the convolution kernel, and the convolution stride is set to 2 to achieve image downsampling. Through strided convolution, image smoothing can be achieved, and some low-frequency rings can be converted into high-frequency rings, achieving better de-ringing effect.

[0048] Step 5: Propose an online dual-domain conversion algorithm. The present invention proposes an online dual-domain conversion algorithm, which can realize the conversion between Cartesian domain and polar coordinate domain while retaining gradient information. Figure 8 shown.

[0049] Cartesian2Polar: Defines the angle θ and radius r in the polar coordinate system. Use torch.meshgrid to combine θ and r into a network that includes the coordinates of each point in polar coordinates. Use the torch.polar function to convert the grid coordinates to polar coordinates, as shown in Equation 3. Finally, use the F.grid_sample method to sample the Cartesian domain image to the polar coordinate domain.

[0050] z=r·cos(θ)+r·sin(θ)·j (4)

[0051] Where r is the radius and θ is the angle (radians).

[0052] Polar2Cartesian: Since the origin of a polar coordinate image is the center of the image, the image is split into a left image and a right image along the central axis, the left image is flipped horizontally, and the two images are spliced ​​vertically. Create a grid, calculate the Cartesian coordinates (X, Y) corresponding to each pixel, and convert the X and Y coordinates to coordinates relative to the center of the image. Calculate the distance r and angle θ from each pixel to the center of the image. Use the F.grid_sample function to sample the spliced ​​image according to the calculated polar coordinate network, obtain the converted image, and adjust the image.

[0053] Step 6: Use the horizontal differential operator to extract the stripe artifacts in the polar coordinate domain and calculate the loss function. The loss function generation process is shown in Fig. 9 After the artifact images at different resolution scales are converted to polar coordinate images through the online coordinate domain conversion algorithm, the horizontal gradient is calculated using the horizontal differential operator. Since the gradient at the stripe artifact is large, the loss function can capture the stripe artifact well. By minimizing the loss function, the stripe artifact is removed in the polar coordinate domain, that is, the ring artifact is removed in the image domain.

[0054] Step 7: Iterate the above process 50 times, the ring artifacts and noise are effectively suppressed, and a ring-free high-quality image is output.

[0055] Example: Figure 3 As shown, the present invention provides an unsupervised inter-layer complementarity enhanced ring artifact removal method, and the ring removal test results are shown in Fig.10 As shown, the specific steps are as follows:

[0056] Step 1: According to the characteristics of the CT scanning system, the ring image is projected into the projection domain to obtain the projection sinusoidal diagram.

[0057] Step 2: Divide the projected sinusoidal image into slice groups, each group includes 5 slices, and each slice group is used as a group of inputs. The input format is [B, S, H, W], where B is the batch size, S is the number of slices, H and W are the height and width of the projected image.

[0058] Step 3: Define the polynomial compensation coefficient equation and initialize the compensation coefficient a j,1 =1, others are 0, calculate the compensated sinusoidal image,

[0059] Step 4: Reconstruct using the FBP reconstruction algorithm and use the TV total variation minimization algorithm to fuse the complementary information between slices for preliminary de-ringing and image denoising.

[0060] Step 5: Use two Gaussian strided convolutions to obtain a multi-scale feature map, and use an online dual-domain conversion algorithm to sample the image into the polar coordinate domain.

[0061] Step 6: Use the horizontal differential operator to calculate the multi-scale feature map function, perform back propagation, and update the compensation coefficient.

[0062] Step 7, iterate the cycle 50 times to obtain a high-quality image without ring artifacts.

[0063] Effect evaluation:

[0064] The present invention discloses an unsupervised ring artifact removal method with enhanced inter-layer complementarity. The test results of the method are as follows: Figure 8 As shown. The figure shows two groups of test images, including a group of model ring artifact data and the de-ringing effect of real ring artifact data. The test was tested in the integrated software. The test process is: import the ring artifact data into the software, then click the "de-ringing button" to get the de-ringed high-quality image, click "Save" to save the de-ringed image in the file, the algorithm runs fast, and it only takes 1.5s to process a slice group. The de-ringing and denoising effects are significant, and the image features can be effectively retained. It shows that this method can quickly and significantly remove ring artifacts, improve image quality, and facilitate subsequent medical and scientific research tasks.

[0065] The above description is only a specific implementation of the present invention, but the protection scope of the present invention is not limited thereto. Any person familiar with the technology can understand and think of any changes or substitutions within the technical scope disclosed by the present invention, which should be included in the scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. An unsupervised ring artifact removal method with inter-layer complementarity enhancement, It is characterized in that The method comprises the following steps: Step 1: According to the characteristics of the CT scanning system, the ring image is projected into the projection domain to obtain the projection sinusoidal diagram. Step 2: Divide the projected sinusoidal image into slice groups, each group includes 5 slices, and each slice group is used as a group of inputs. The input format is [B, S, H, W], where B is the batch size, S is the number of slices, H and W are the height and width of the projected image. Step 3: Define the polynomial compensation coefficient equation and initialize the compensation coefficient a j,1 =1, the others are 0, calculate the compensated sinusoidal image, step 4, reconstruct using the FBP reconstruction algorithm, use the TV total variation minimization algorithm to fuse the complementary information between slices for preliminary de-ringing and image denoising, Step 5: Use two Gaussian strided convolutions to obtain a multi-scale feature map, and use an online dual-domain conversion algorithm to sample the image into the polar coordinate domain. Step 6: Use the horizontal differential operator to calculate the multi-scale feature map function, perform back propagation, and update the compensation coefficient. Step 7, iterate the cycle 50 times to obtain a high-quality image without ring artifacts.

2. The unsupervised ring artifact removal method with inter-layer complementarity enhancement according to claim 1, It is characterized in that The projection mentioned in step 1 is performed according to the characteristics of the CT scanning system, and the system characteristics include the reconstructed image size, the number of detector pixels, the detector pixel size, the detector scanning angle, and the distance between the detector and the object.

3. The unsupervised ring artifact removal method with inter-layer complementarity enhancement according to claim 1, It is characterized in that In step 2, the slice group input format, that is, the 2.5D input format, is based on the principle that the detector has randomness in the head and foot directions, resulting in the complementarity of the annular artifact information of each slice. Using 2.5D data as input can retain the spatial redundant information between slices.

4. The unsupervised ring artifact removal method with inter-layer complementarity enhancement according to claim 1, It is characterized in that The polynomial compensation coefficient equation in step 3, the main reason for the ring artifact is the non-uniform response of the detector. By optimizing the detector response, the ring artifact can be effectively removed. The specific method is: define {a j,0 ,a j,1 ,...,a j,i }Detector j-th pixel value S j The compensation coefficient of 0~i power, compensation equation Optimize j,0 ,a j,1 ,...,a j,i } Positive detector response, achieving the purpose of removing ring artifacts.

5. The unsupervised ring artifact removal method with inter-layer complementarity enhancement according to claim 1, It is characterized in that In step 4, the TV total variation algorithm extracts inter-slice complementary information and denoises. The details are as follows: the total variation of the reconstructed image in the Z, X, and Y directions is calculated using the following formula, and a fidelity term is added to prevent large-area distortion of the image during the TV minimization process. After each iteration, the image is added to the normalized total variation to smooth the noise and eliminate some ring artifacts. img corr =img+Norm(J) Among them, the first is the total variation of the image. Limiting the total variation can suppress noise. The second term is the fidelity term, which prevents the denoised image from being too different from the original image. λ is the relaxation factor. Since K is small, the image distortion phenomenon is not obvious, so it is set to λ=0 in this algorithm.

6. The unsupervised inter-layer complementarity enhanced ring artifact removal method according to claim 1, It is characterized in that The online dual-domain conversion algorithm in step 5 is as follows: Cartesian2polar: defines the polar coordinate domain radius r and angle θ, generates a polar coordinate grid, and uses the F.grid_sample method to resample the image to the polar coordinate system. Polar2cartesian: Defines Cartesian coordinates X and Y, generates a Cartesian coordinate network, splits and reassembles the image along the central axis, calculates the Cartesian coordinates (X, Y) corresponding to each pixel, converts the X and Y coordinates to coordinates relative to the center of the image, converts the coordinates to polar coordinates, and uses F.grid_sample to resample the image to the Cartesian domain.

7. The unsupervised ring artifact removal method with inter-layer complementarity enhancement according to claim 1, It is characterized in that The Gaussian strided convolution mentioned in step 5 downsamples the image using a Gaussian kernel as the convolution kernel. The Gaussian kernel can smooth the image and reduce some of the image noise. Using strided convolution, some low-frequency ring artifacts are sampled as high-frequency ring artifacts, thereby enhancing the artifact removal effect.

8. The unsupervised inter-layer complementarity enhanced ring artifact removal method according to claim 1, It is characterized in that The horizontal differential operator mentioned in step 6 calculates the gradient. The ring artifact appears as vertical stripes in the polar coordinate domain. The horizontal differential can effectively extract the artifact gradient under the multi-scale feature map, and the compensation coefficient is optimized through back propagation.

Images