Method for using weighted multi-scale local anisotropic total variation based on prior guidance for CT image Missing Wedging artifact correction

Abstract

The invention relates to a multi-scale local anisotropic total variation (PG-WM-LATV) method based on prior guidance, and the method comprises the following steps: constructing a composite objective function in combination with a prior image, including a weighted data fidelity term and a multi-directional regularization term, the method comprises the following steps of: firstly, performing multi-scale post-processing on a target function, continuously updating dual variables and original variables of the target function through a Chambolle-Pock (CP) algorithm optimization framework, then performing multi-scale post-processing, decomposing a to-be-processed image into a plurality of scales by a Gaussian pyramid in the multi-scale post-processing, and realizing frequency domain compensation, local projection domain restoration and direction regularization at each scale. Compared with an experiment of an existing reconstruction method, the method has the advantages that region information missing due to Missing Wedging can be effectively recovered, and the image structure fidelity and the detail recovery effect are remarkably improved. And an effective solution is provided for the Missing Wedging problem in the field of CT imaging.

Description

Technical Field

[0001] The present invention relates to a MissingWedge artifact processing method based on a priori-guided weighted multi-scale local anisotropy total variation, and belongs to the field of medical CT, industrial CT and material CT imaging. Background Art

[0002] Computed tomography (CT), which reconstructs the internal three-dimensional structure of an object from projection data, has become a core technology in various fields, including material microscopic characterization, modern medical diagnosis, and industrial non-destructive testing. Its core advantage lies in its ability to provide three-dimensional information about the target's internal structure. However, in actual imaging, due to hardware conditions, object imaging state, radiation damage, and dose limitations, tilt angle sampling is incomplete. According to the Fourier slice theorem, unsampled high-tilt angle regions form wedge-shaped gaps in Fourier space, known as "missing wedges." This leads to anisotropic resolution and wedge-shaped artifacts, seriously compromising the accuracy of image reconstruction and the reliability of subsequent analysis.

[0003] To address the Missing Wedge problem, various data analysis algorithms have been proposed for the missing information. (1) Traditional reconstruction algorithms. The Filtered Back Projection (FBP) algorithm is a widely used tomography method, but it is extremely sensitive to Missing Wedge and can easily cause streak artifacts and axial stretching. The Weighted Back Projection (WBP) algorithm can compensate for some missing frequencies by introducing an angle-dependent weight function, but there are obvious wedge artifacts when the tilt angle is <±60°. The Direct Fourier Method (DFM) directly reconstructs in Fourier space through polar coordinate-Cartesian interpolation, but the interpolation error is amplified in the Missing Wedge area, resulting in ring artifacts. (2) Iterative reconstruction algorithm. Classic iterative reconstruction algorithms include Simultaneous Iterative Reconstruction Technique (SIRT), Simultaneous Algebraic Reconstruction Technique (SART), and Discrete Algebraic Reconstruction Technique (DART). These algorithms are based on the iterative optimization framework of linear models and gradually correct the image estimate by minimizing the difference between the original projection and the projection of the reconstructed volume. However, they are sensitive to noise and cannot recover the real information within the Missing Wedge. (3) Regularization (Total Variation, TV) algorithms based on compressed sensing. For example, the Total Variation Minimization (TVM) method solves the incomplete projection problem by minimizing the image gradient lp norm. However, TVM will cause the real edge details to be over-smoothed. To overcome this limitation, many advanced TV reconstruction methods have been proposed, such as the reweighted anisotropic total variation (RwATV) algorithm, the adaptive weighted TV (AwTV) algorithm, the TV algorithm based on the exponentially directional weighted function (EDWF), the directional TV (DTV) algorithm, and the local anisotropic total variation (LATV) algorithm.These algorithms can achieve better reconstruction accuracy and faster convergence speed under limited angles, but a single-scale framework has difficulty in synergistically recovering low-frequency global structures and high-frequency local details.

[0004] In recent years, the explosive development of machine learning methods such as deep neural networks has made significant contributions to high-performance information recovery for the Missing Wedge problem. On the one hand, neural networks can further optimize traditional reconstruction algorithms. For example, Deep FBP and PCNet maintain high computational efficiency while significantly reducing sampling costs. On the other hand, neural networks can overcome the limitations of traditional algorithms, such as the UsiNet (Unsupervised Sinogram Inpainting for Nanoparticle Electron Tomography) method. Although deep learning is fast at the inference stage, it relies on large amounts of high-quality paired data, and hardware resource limitations can hinder its deployment.

[0005] Therefore, to address the shortcomings of the existing technology, the present invention proposes a prior-guided weighted multiscale local anisotropic total variation (PG-WM-LATV) method. Summary of the Invention

[0006] The present invention provides a PG-WM-LATV method for Missing Wedge artifact correction in CT images. This method comprises the following steps: constructing a composite objective function based on a prior image, comprising a weighted data fidelity term and a multi-directional regularization term; applying the objective function to the Chambolle-Pock (CP) algorithm optimization framework to continuously update the dual and primal variables; and then performing multi-scale post-processing. This multi-scale post-processing decomposes the image to be processed into multiple scales using a Gaussian pyramid, implementing frequency domain compensation, local projection domain restoration, and directional regularization at each scale. This method can effectively restore regional information lost due to Missing Wedges and significantly improve image structural fidelity and detail recovery. The present invention comprises the following steps:

[0007] Step 1: Reconstruct a high-quality image using projection data from a full or nearly full-angle scan, using this as the prior image. Based on the Parker function's constraint on weights, a weighting function is designed such that the sum of the two weights corresponding to the same line integral is 1, and the weights in areas with missing data are 0. Furthermore, the weights themselves and their gradients must be continuous across the entire 360° range, achieving a smooth transition at the missing data boundaries. This avoids spectral discontinuities caused by binary weighting and effectively suppresses artifacts.

[0008] Step 2: The Chambolle-Pock (CP) algorithm optimization framework is used to facilitate subsequent experimental comparisons, and is called PG-W-LATV. The CP algorithm can be used to perform CT image reconstruction on a specific convex optimization problem, converting the objective function into a saddle point problem. The following steps are followed to perform alternating updates:

[0009] Step 2.1: Update the data fidelity term dual variable w. Update the closed-form solution obtained based on the proximal operator and the weighted residual.

[0010] Step 2.2: Update the dual variables of the multi-directional regularization term . Update according to the projection operator projected to the l1 sphere.

[0011] Step 2.3: Update the original variable u. Update the original variable according to steps 2.1 and 2.2.

[0012] Step 2.4: Acceleration step, introduce extrapolation to improve the convergence speed.

[0013] Repeat steps 2.1 to 2.4 until convergence is achieved.

[0014] Step 3: Multi-scale post-processing. The multi-scale post-processing strategy is based on LATV optimization and divides the regularization process into multi-scale decomposition and single-scale iteration. Through multi-scale decomposition, the low-frequency, medium-frequency, and high-frequency components of the image are gradually restored. Specifically, the image is first decomposed into multi-scale residuals. Then, frequency domain compensation, local projection domain repair, and direction regularization are performed within each scale. Finally, the low-resolution result is upsampled and fused with the high-frequency residual. The following steps are followed to update alternately:

[0015] Step 3.1: Multi-scale decomposition: Construct a Gaussian pyramid and perform weighted fusion of the preliminary reconstruction result and the original pyramid residual.

[0016] Step 3.2: Single-scale iteration. This step is an inner loop consisting of three parts: frequency domain compensation, local projection domain restoration, and directional TV regularization. Frequency domain compensation extracts local blocks from missing points, calculates the gradient matrix, and then weightedly fuses the frequency domains of adjacent angles. Local projection domain restoration adjusts the projection values ​​along the gradient direction. Directional TV regularization produces the updated image after directional TV regularization.

[0017] The beneficial effects of the present invention are as follows: (1) The present invention simultaneously introduces a weighted data fidelity term and a multi-directional regularization term into the objective function. This effectively suppresses the Fourier spectrum discontinuity caused by data loss and enhances the image edge recovery capability. (2) Based on the CP algorithm, the present invention designs a dual variable update formula and a primal variable update strategy, and ensures step size update by calculating the operator norm, thereby achieving efficient and stable convergence in large-scale CT reconstruction. (3) The multi-scale post-processing strategy of the present invention has significant advantages in image resolution and detail recovery, enabling the PG-WM-LATV algorithm to still exhibit superior performance under conditions of incomplete data and noise interference. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the present invention is further described in detail below with reference to the accompanying drawings and specific implementation methods:

[0019] Figure 1 : Flowchart of the method of the present invention.

[0020] Figure 2 : Comparison of reconstruction results of five methods of Shepp-Logan model under different projection angles.

[0021] Figure 3 : Figure 2 Comparison of SSIM, PSNR, MSE and IoU values ​​of different reconstruction algorithms under corresponding projection angles.

[0022] Figure 4 :The Shepp-Logan model is reconstructed using different algorithms at signal-to-noise ratios of 20dB, 25dB, and 30dB. The projection angle ranges from 0° to 140°.

[0023] Figure 5 : Schematic diagram of the absorption contrast imaging experimental setup.

[0024] Figure 6 : CT scanning experimental data of the head of dried mealworm samples.

[0025] Figure 7 Comparison of reconstruction results of mealworm head slices using different algorithms. Projection angles range from 0° to 120°. DETAILED DESCRIPTION

[0026] In order to enable those skilled in the art to better understand the present application, the present invention will be further described below with reference to the accompanying drawings and specific implementation methods.

[0027] like Figure 1 As shown, the PG-WM-LATV method of the present invention comprises the following steps:

[0028] Step 1: First, use the projection data of the full angle or almost full angle scan to reconstruct a high-quality image. Based on the constraint principle of Parker function on weights, design a weighting function so that the sum of the two weights corresponding to the same line integral is 1, and the weight of the area with missing data is 0. The weight itself and the gradient of the weight must be continuous within the entire 360° range, and a smooth transition is achieved at the missing boundary to avoid spectrum discontinuity caused by binary weighting and effectively suppress artifacts. Use it as a priori image: Let the image to be reconstructed be u, and observe the projection data , construct a composite objective function : (1)

[0029] Step 2: The Chambolle-Pock (CP) algorithm optimization framework is used to facilitate subsequent experimental comparisons and is called PG-W-LATV. The CP algorithm solves the following saddle point optimization problem by alternating between updating the primal and dual variables. The CP algorithm can be used to perform CT image reconstruction on specific convex optimization problems. By converting the objective function into a saddle point problem, we obtain: (2) Among them, w is the dual variable of the weighted data fidelity term, and qk is the dual variable of the multi-directional regularization term.

[0030] Follow these steps to update alternately:

[0031] Step 2.1: Update the data fidelity term dual variable w. Update the closed-form solution obtained based on the proximal operator and weighted residual: (3)

[0032] Step 2.2: Update the dual variables of the multi-directional regularization term . Update according to the projection operator projected to the l1 sphere: (4)

[0033] Step 2.3: Update the original variable u. Update the original variable according to steps 2.1 and 2.2.

[0034] Step 2.4: Acceleration step, introduce extrapolation to improve the convergence speed.

[0035] Step 3: Multi-scale post-processing. The multi-scale post-processing strategy is based on LATV optimization and divides the regularization process into multi-scale decomposition and single-scale iteration. Through multi-scale decomposition, the low-frequency, medium-frequency, and high-frequency components of the image are gradually restored. Specifically, the image is first decomposed into multi-scale residuals. Then, frequency domain compensation, local projection domain repair, and direction regularization are performed within each scale. Finally, the low-resolution result is upsampled and fused with the high-frequency residual. The following steps are followed to update alternately:

[0036] Part 1: Frequency domain compensation extracts local blocks for missing points (r, t) and calculates the gradient matrix, and weighted fusion of the frequency domain F (k+1) of adjacent angles is: (5) Where h(k) is the Hanning window function.

[0037] Part 2: Local projection domain repair, adjusting the projection value along the gradient direction: (6) Among them, the tanh function can control the update amplitude.

[0038] Part 3: Directional TV regularization, u(n) is the reconstructed image after the nth iteration, u(n+1) is the image updated after directional TV regularization, then: (7)

[0039] Repeat these three parts until all scales are completed.

[0040] We verify the effectiveness of the implemented method through experimental data of the Shepp-Logan model and head CT scans of dried mealworm samples.

[0041] like Figure 2 As shown in the figure, the Shepp-Logan model is used as the test object to verify the robustness of this method. Figure 2 Reconstruction results of five methods of Shepp-Logan model at different projection angles. (a)-(e) are projections from 0° to 120°, (f)-(j) are projections from 0° to 140°, (k)-(o) are projections from 0° to 160°, (a), (f) and (k) are TV methods, (b), (g) and (l) are RwATV methods, (e), (h) and (m) are LATV methods, (d), (i) and (n) are PG-W-LATV methods, (e), (j) and (o) are PG-WM-LATV methods. Figure 2 It can be seen that at different projection angles, the wedge-shaped artifacts in the reconstructed image using the PG-WM-LATV method almost disappear, and the boundary details are most significantly restored.

[0042] Figure 3 for Figure 2 Comparison of SSIM, PSNR, MSE, and IoU values ​​for different reconstruction methods at corresponding projection angles. (a)-(c) Scatter plots of the SSIM and PSNR for the reconstructions from 0° to 120°, 0° to 140°, and 0° to 160°, respectively, with the horizontal axis representing SSIM and the vertical axis representing PSNR. (d)-(f) Scatter plots of the MSE and IoU for the reconstructions from 0° to 120°, 0° to 140°, and 0° to 160°, respectively, with the horizontal axis representing MSE and the vertical axis representing IoU. Experimental data shows that the PG-WM-LATV method achieves significantly higher SSIM and PSNR values ​​at different projection angles, with SSIM improvements of approximately 0.35-0.37 and PSNR improvements of approximately 7.1dB-8.3dB compared to the TV method. Moreover, this method has the smallest MSE value and the highest IoU value, which shows that the PG-W-LATV method has the best performance and can reconstruct edge and region information more accurately and restore the original structure more accurately.

[0043] Figure 4 It can better verify the anti-noise ability of the proposed method. Figure 4 Different methods were used to reconstruct the Shepp-Logan model at signal-to-noise ratios of 20dB, 25dB, and 30dB. (a)-(c) 20dB, (d)-(f) 25dB, and (g)-(i) 30dB. (a), (d), and (g) LATV methods, (b), (e), and (h) PG-W-LATV methods, and (c), (f), and (i) PG-WM-LATV methods. The Shepp-Logan model projection range is 0° to 140°. As can be seen from the images, even in the presence of noise, the PG-WM-LATV method reconstructs images that are relatively smooth and restores both structural and edge information, demonstrating its strong noise immunity.

[0044] Figure 5 Schematic diagram of the absorption contrast imaging experimental device. In order to more realistically verify the practicality and effectiveness of this method, absorption contrast imaging was used to obtain a set of head CT scan experimental data of dried mealworm samples.

[0045] Figure 6 The following is the experimental data of CT scanning of the head of a dried mealworm sample. The experiment obtained 600 equal-angle projections from the dried mealworm sample. The rotation angle was 0° to 179.9°. The X-ray photon energy was 12keV. The exposure time of each projection was 0.055 s. The bright field background image and dark field background image were collected at the same time for preprocessing the obtained projections. The preprocessed data is as follows: Figure 6 shown. Figure 6(a) is a schematic diagram of the three-dimensional projection, (b)-(i) are the 50th, 125th, 200th, 275th, 350th, 425th, 500th and 575th projection images after preprocessing, with a pixel size of 573×1024. Figure 6 It can be seen that this data has high resolution and rich tissue structure details, which can effectively verify the application of the PG-WM-LATV method in the field of biological imaging.

[0046] Figure 7 Comparison of the reconstruction effects of different methods on the slices of the mealworm head. The 600 pre-processed projections were reconstructed to obtain 573 slices, such as Figure 7 (a) Original full-angle CT scan image, 950×950 pixels; (b) LATV reconstruction. (c) PG-W-LATV reconstruction. (d) PG-WM-LATV reconstruction. Projection angles range from 0° to 120°. The PG-WM-LATV reconstruction shows the most significant restoration of overall image structure, edge details, and texture information. Wedge-shaped artifacts are largely eliminated, and image details are still restored even in low-contrast areas.

[0047] The above describes the basic principles, main features and advantages of the present invention. It should be understood that those skilled in the art can improve the present invention based on the above description, and all such improvements and modifications should be included in the scope of the claims of the present invention.

Claims

1. A method for Missing Wedge artifact correction in CT images based on prior-guided weighted multi-scale local anisotropy total variation, characterized by The following steps are involved: Step 1: Use projection data from full-angle or nearly full-angle scans to reconstruct a high-quality image and use it as a priori image. Based on the Parker function's constraint on weights, a weighting function is designed so that the sum of the two weights corresponding to the same line integral is 1, and the weight of the area with missing data is 0. Furthermore, the weight itself and its gradient must be continuous throughout the entire 360° range, achieving a smooth transition at the missing data boundary to avoid spectral discontinuity caused by binary weighting and effectively suppress artifacts. Step 2: The Chambolle-Pock (CP) algorithm optimization framework is used to facilitate subsequent experimental comparisons, and is called PG-W-LATV. The CP algorithm can be used to perform CT image reconstruction on a specific convex optimization problem, converting the objective function into a saddle point problem. The following steps are followed to perform alternating updates: Step 2.1: Update the data fidelity term dual variable w. Update the closed-form solution obtained based on the proximal operator and the weighted residual. Step 2.2: Update the dual variables of the multi-directional regularization term . Update according to the projection operator projected to the l1 sphere; Step 2.3: Update the original variable u. Update the original variable according to steps 2.1 and 2.2; Step 2.4: Acceleration step, introducing extrapolation to improve the convergence speed; Repeat steps 2.1 to 2.4 until convergence is achieved; Step 3: Multi-scale post-processing. The multi-scale post-processing strategy is based on LATV optimization and divides the regularization process into multi-scale decomposition and single-scale iteration. Through multi-scale decomposition, the low-frequency, medium-frequency, and high-frequency components of the image are gradually restored. Specifically, the image is first decomposed into multi-scale residuals. Then, frequency domain compensation, local projection domain repair, and direction regularization are performed within each scale. Finally, the low-resolution result is upsampled and fused with the high-frequency residual. The following steps are followed to update alternately: Step 3.1: Multi-scale decomposition. Construct a Gaussian pyramid and perform weighted fusion of the initial reconstruction result and the original pyramid residual. Step 3.2: Single-scale iteration. This step is an inner loop consisting of three parts: frequency domain compensation, local projection domain restoration, and directional TV regularization. Frequency domain compensation extracts local blocks from missing points, calculates the gradient matrix, and then weightedly fuses the frequency domains of adjacent angles. Local projection domain restoration adjusts the projection values ​​along the gradient direction. Directional TV regularization produces the updated image after directional TV regularization.

2. The method according to claim 1, characterized in that Composite objective function The formula is:

3. The method according to claim 1, characterized in that The CP algorithm can be used to reconstruct CT images for specific convex optimization problems. The formula for converting the objective function into a saddle point problem is: Among them, w is the dual variable of the weighted data fidelity term, and qk is the dual variable of the multi-directional regularization term.

4. The method according to claim 1, wherein The formula for updating the data fidelity term dual variable w is:

5. The method according to claim 1, wherein Update the dual variables of the multi-directional regularization term . Update according to the projection operator projected to the l1 sphere:

6. The method according to claim 1, characterized in that Frequency domain compensation extracts local blocks for missing points (r, t) and calculates the gradient matrix, and weighted fusion of the frequency domain F (k+1) of adjacent angles is: Where h(k) is the Hanning window function.

7. The method according to claim 1, characterized in that Repair the local projection domain and adjust the projection value along the gradient direction: Among them, the tanh function can control the update amplitude.

8. The method according to claim 1, characterized in that Directional TV regularization, u(n) is the reconstructed image after the nth iteration, u(n+1) is the image updated after directional TV regularization, then:

Images