Sparse angle CT image reconstruction method based on compressed sensing theory

Abstract

The invention discloses a sparse angle CT image reconstruction method based on a compressed sensing theory. According to the scheme, an FBP initial reconstruction image is input into an ART reconstruction module; storing the current image; traversing all angle numbers; extracting the current angle, calculating the forward projection of the current angle, and calculating and storing the projection error of the current angle; the angle internal circulation is ended; performing back projection on the projection error at the current angle and updating the current image; fISTA denoising is executed, and the current image is updated; and judging whether convergence occurs or not, if not, continuing the ART iteration module, and ending the ART iteration module. According to the scheme, a novel compressed sensing algorithm, an HWDSC algorithm and an AWL1R algorithm are designed, so that artifacts and noise are remarkably reduced under the sparse sampling condition, and the radiation dose of a patient is reduced.

Description

Technical Field

[0001] The present invention relates to the field of medical image processing, and particularly to a sparse-angle CT image reconstruction method based on the theory of compressive sensing. Background Art

[0002] Since the advent of CT technology in the 1970s, it has become the core imaging means in medical diagnosis and industrial inspection. Its basic principle is to penetrate the target object with X-rays, collect projection data from multiple angles, and use mathematical algorithms to reconstruct the internal structure image of the target. However, in order to ensure high resolution and low noise, traditional CT scans usually need to collect a large amount of projection data. For example, more than 180 projections are made in the range of 0° to 179°. This high sampling rate exposes patients to a higher radiation dose. Research shows that excessive radiation may increase health risks such as cancer and leukemia, especially having a more significant impact on sensitive populations such as children and pregnant women. Therefore, how to maintain image quality under the condition of reducing the number of projections (i.e., sparse-angle sampling) has become a key research direction in the field of CT technology.

[0003] In traditional reconstruction methods, Filtered Back Projection (FBP), as a classic analytical algorithm, generates an image by applying a filter (such as the Ram-Lak filter) to the projection data and back-projecting it. It has the advantages of high computational efficiency and good results under fully sampled conditions. However, when the number of projections is reduced to less than 60 times, serious strip artifacts and noise will appear in the FBP reconstructed image, resulting in blurred edges and lost details, making it difficult to meet the accuracy requirements of clinical diagnosis. Algebraic Reconstruction Technique (ART) can alleviate the artifact problem caused by sparse sampling to a certain extent by iteratively optimizing the projection consistency, but it is sensitive to the initial value and noise, has a slow convergence speed, and the reconstruction quality is still not ideal under extremely sparse conditions (such as when the number of projections is less than 30 times).

[0004] In recent years, the CS theory has provided a new theoretical basis for sparse-angle CT reconstruction. CS points out that if a signal is sparse in a certain transform domain, the original signal can be accurately reconstructed from undersampled data. This theory has been widely applied to low-dose CT imaging, and common methods include Total Variation (TV) regularization and Wavelet Transform. TV regularization effectively suppresses noise and artifacts by minimizing the image gradient, but it is prone to over-smoothing the image and losing edges and fine structures under sparse sampling conditions; the wavelet transform is good at capturing local features of images, but it has limited ability to represent global structures and may introduce blocky artifacts. Existing CS methods mostly rely on a single sparse transform domain and cannot simultaneously take into account the local details and global consistency of images. Especially under complex structures (such as the junction of bone and soft tissue) or extremely sparse sampling conditions, the reconstruction effect is not good.

[0005] Current research shows that the limitations of a single sparse transform can be compensated by multi-domain sparse representation. For example, combining the wavelet transform and the Discrete Cosine Transform (DCT) can improve sparsity. However, how to effectively fuse multi-domain features and design an adaptive optimization strategy to adapt to different sampling conditions, noise levels, and image complexities remains an unsolved technical problem. In addition, the regularization parameters of traditional methods are usually fixed values and lack the ability of dynamic adjustment, resulting in insufficient robustness of the algorithm in different scenarios. Summary of the Invention

[0006] In view of the above deficiencies in the prior art, the present invention provides a sparse-angle CT image reconstruction method based on the compressed sensing theory.

[0007] To achieve the above invention purpose, the technical solution adopted by the present invention is as follows: A sparse-angle CT image reconstruction method based on the compressed sensing theory, comprising the following steps: S1. Load projection data and initialize parameter settings, where the parameters include the number of projections, the number of iterations, and the tolerance; S2. Apply the FBP method to the projection data and output an initial reconstructed image; S3. Input the obtained initial reconstructed image into the ART reconstruction module and perform ART module iterations; S4. Perform external optimization on the image output after ART module iterations, and update the projection image after the optimization ends.

[0008] Further, the S3 specifically includes the following steps: S31. Input the obtained initial reconstructed image into the ART reconstruction module and save the current image; S32, traverse all angle numbers, extract the current angle and calculate the forward projection of the current angle, calculate and save the current angle projection error, and end the angle inner loop; S33, performing backprojection on the projection error at the current angle and updating the current image, performing FISTA denoising and updating the current image; S34: Determine whether the updated image has converged. If not, continue to iterate the ART module until it converges.

[0009] Furthermore, the specific formula of the ART module in S34 is:

[0010] Where x is the image to be reconstructed, a i is the i-th row of the projection matrix A, y i is the projection data, λ is the relaxation parameter, k is the number of iterations, T is the matrix transpose, Calculate the norm.

[0011] Furthermore, the external optimization in S4 adopts an optimization algorithm based on hybrid wavelet and DCT coefficient constraint or an optimization algorithm based on adaptive wavelet transform and L1 regularization.

[0012] Furthermore, the optimization algorithm based on hybrid wavelet and DCT coefficient constraints takes the reconstructed image of filtered back projection as the initial value, and its optimization objective function is expressed as:

[0013] Where, x is the wavelet transform operator, λ1 and λ2 are constants, x is the image to be reconstructed, DCT is the discrete cosine transform operator, is the norm calculation, y is the measurement data, and A is the projection matrix.

[0014] Furthermore, the optimization algorithm based on adaptive wavelet transform and L1 regularization takes the reconstructed image of filtered back projection as the initial value, and its optimization objective function is expressed as:

[0015] Where x is the image to be reconstructed, μ is the L1 regularization parameter in the image domain, is the norm calculation, y is the measurement data, A is the projection matrix; λ is the wavelet domain regularization parameter, and the adaptive threshold adjustment method is: ; Where a is a constant, n is the current number of projections, N is the maximum number of projections, i and I are the current and total number of iterations, respectively.

[0016] The present invention has the following beneficial effects: 1. Compared with the traditional FBP algorithm, this method significantly reduces artifacts and noise under sparse sampling conditions, improves objective metrics such as Peak Signal-to-Noise Ratio (PSNR), Structural Similarity Index (SSIM), and Mean Squared Error (MSE), and simultaneously enhances the subjective visual effect.

[0017] 2. This method enhances the algorithm's robustness: through multi-domain sparse representation and adaptive regularization, it ensures the stability and adaptability of the algorithm under different projection numbers (20 - 180 times), noise levels (0% - 10%), and image complexities (such as bone and soft tissue mixed regions).

[0018] 3. This method supports low-dose CT imaging: it reduces 50% - 90% of the projection data, lowers the patient's radiation dose, and simultaneously ensures that the reconstructed images meet the high-precision requirements of clinical diagnosis and industrial inspection.

[0019] 4. Compared with deep learning algorithms with relatively large computational requirements, this method has the advantage of optimized computational efficiency. Under the premise of ensuring image quality, it shortens the reconstruction time (targeting 5 - 10 seconds) to meet the real-time imaging requirements. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 It is a schematic diagram of the algorithm framework of the sparse-angle CT reconstruction method based on the compressed sensing theory proposed by the present invention.

[0021] Figure 2 It is a reconstructed image of the FBP sampling projection number in the embodiment of the present invention.

[0022] Figure 3 It is a reconstructed image of the ART sampling projection number in the embodiment of the present invention.

[0023] Figure 4 It is a reconstructed image of the HWDSC sampling projection number in the embodiment of the present invention.

[0024] Figure 5 It is a reconstructed image of the AWL1R sampling projection number in the embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0025] The following describes the specific embodiments of the present invention to facilitate those skilled in the art to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions created using the concept of the present invention are within the scope of protection.

[0026] A sparse angle CT image reconstruction method based on compressed sensing theory, such as Figure 1 As shown, the following steps are included: S1, load the projection data and initialize the parameter settings, where the parameters include the number of projections, the number of iterations and the tolerance; To validate the effectiveness of the proposed algorithm, experiments were conducted using a 256×256 pixel Shepp-Logan skull phantom. Projection data was generated using a Radon transform, with the number of sampled projections set to 180, 90, 60, 30, and 20 within the 0° to 179° range to simulate a sparse sampling scenario. This standard test image, widely used in CT reconstruction research, features piecewise smooth regions and sharp edges, making it suitable for evaluating the performance of reconstruction algorithms. To quantify the reconstruction performance, three objective metrics were designed and calculated: mean squared error (MSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM). MSE measures the pixel-level difference between the reconstructed image and the original image, with smaller values indicating better quality. PSNR reflects the signal-to-noise ratio (SNR) of the image, measured in decibels (dB), with higher values indicating better quality. SSIM assesses the structural similarity of images, with values ranging from 0 to 1, with values closer to 1 indicating higher similarity.

[0027] S2, applying the FBP method to the projection data and outputting an initial reconstructed image; S3, inputting the obtained initial reconstructed image into the ART reconstruction module, and performing ART module iteration; S31, inputting the obtained initial reconstructed image into the ART reconstruction module and saving the current image; S32, traverse all angle numbers, extract the current angle and calculate the forward projection of the current angle, calculate and save the current angle projection error, and end the angle inner loop; S33, performing backprojection on the projection error at the current angle and updating the current image, performing FISTA denoising and updating the current image; ART iterative module: input the FBP initial reconstructed image into the ART reconstruction module; save the current image; then traverse all angles; extract the current angle, calculate the current angle forward projection, calculate and save the current angle projection error; end the angle inner loop; perform back projection on the projection error at the current angle and update the current image; perform FISTA (internally once) denoising and update the current image; determine whether it converges. If not, continue with the ART iterative module and end the ART iterative module.

[0028] Among them, the ART iterative update formula is: (1) Where x∈R n is the image to be reconstructed, a iis the i-th row of the projection matrix A, and y i is the projection data, and λ ∈ (0, 2) is the relaxation parameter (a typical value is 0.1).

[0029] S34. Determine whether the updated image converges. If it does not converge, continue to execute the ART module iteration until it converges.

[0030] S4. Externally optimize the image output after the ART module iteration. After the optimization is completed, update the projection image.

[0031] External optimization module: Externally apply the FISTA algorithm to the image output by the ART iteration module, and execute one of the following two optimization algorithms, including the optimization algorithm based on the hybrid wavelet and DCT coefficient constraint (HWDSC) or the optimization algorithm based on the adaptive wavelet transform and L1 regularization (AWL1R). After the optimization is completed, update the image.

[0032] Among them, the principle of the HWDSC (hybrid wavelet + DCT sparse constraint) algorithm is to combine the multi-domain sparse constraints of wavelet transform (WaveletTransform) and discrete cosine transform (Discrete Cosine Transform, DCT), and make full use of the local sparsity of the image in the wavelet domain and the global frequency domain sparsity in the DCT domain. The multi-domain sparse constraints effectively balance the local details (such as edges and textures) and global smoothness (such as uniform tissue regions), significantly reducing artifacts and improving image consistency. The optimization objective function of the HWDSC algorithm is:[[]] (2) Among them, Ψ x is the wavelet transform operator (using the 4-level Daubechies wavelet), λ1 = 0.0004, and λ2 = 0.00005.

[0033] The principle of the AWL1R (adaptive wavelet + L1 regularization) algorithm is as follows: Introduce the adaptive wavelet threshold and L1 regularization in the image domain, and dynamically adjust the regularization parameter to adapt to different sampling conditions and iteration stages. The adaptive regularization mechanism enhances the robustness of the algorithm under extreme sparse conditions (such as 20 projections), and the L1 regularization in the image domain further suppresses noise and improves the stability of smooth regions. The optimization objective function of the AWL1R algorithm is:[[]] (3) Among them, μ = 0.00005 is the L1 regularization parameter in the image domain. λ is the regularization parameter in the wavelet domain, and its adaptive threshold adjustment method is:[[]] (4) where a = 0.01, n is the current number of projections, N = 180 is the maximum number of projections, i and I = 35 are the current and total number of iterations respectively, and 10 -6 Avoid division by zero.

[0034] Image quality assessment module: objective metrics such as Peak Signal-to-Noise Ratio (PSNR), Structural Similarity Index (SSIM), and Mean Squared Error (MSE), output the metrics and the image, and end the reconstruction algorithm.

[0035] Parameter settings in the experiment: For the HWDSC algorithm, the FISTA regularization parameter inside ART is 0.00015, the external FISTA regularization parameter (in the wavelet domain) is 0.0004, the sparse weight in the DCT domain is 0.00005, and the number of external FISTA iterations is 30 times.

[0036] For the AWL1R algorithm, the parameters are: the FISTA regularization parameter inside ART is 0.00015, the external FISTA regularization parameter (based on the wavelet domain) is 0.0004, the L1 regularization weight is 0.0001, and the number of external FISTA iterations is 35 times.

[0037] The experiment compares the reconstruction effects of FBP (using Hann filtering), ART (50 iterations), and the two CS methods proposed in the present invention. The results are shown in Table 1. Under extremely sparse sampling conditions, when the number of projections is 20, the PSNR of the HWDSC optimization method is 26.59 dB and the SSIM is 0.9191; the PSNR of the AWL1R optimization method is 26.38 dB and the SSIM is 0.9166, both significantly superior to FBP and ART. When the number of sampling projections is 20, the images show that the CS methods can effectively reduce artifacts and recover details.

[0038] Table 1 Comparison of reconstruction quality

[0039] The beneficial effects of the present invention are mainly reflected in the following aspects. Through multi-domain sparse optimization and adaptive regularization techniques, combined with the simulation experimental data as shown in Table 1, the results are as Figures 2 - 5 shown, significantly improving the performance of sparse-angle CT reconstruction. The following elaborates on each advantage in detail: (1) Significantly improve image quality: The HWDSC (Hybrid Wavelet + DCT Sparse Constraint) and AWL1R (Adaptive Wavelet + L1 Regularization) algorithms proposed by the present invention significantly improve the image quality under sparse angular conditions, specifically manifested as excellent performance in objective indicators such as Peak Signal-to-Noise Ratio (PSNR), Structural Similarity (SSIM), and Mean Squared Error (MSE).

[0040] Under fully sampled conditions (180 projections): According to Table 1, the PSNR of the HWDSC algorithm reaches 30.31 dB and the SSIM is 0.9668; the PSNR of the AWL1R algorithm is even higher at 31.09 dB and the SSIM is 0.9754, far exceeding the traditional Filtered Back Projection (FBP, PSNR = 24.59 dB, SSIM = 0.8807) and Algebraic Reconstruction Technique (ART, PSNR = 23.67 dB, SSIM = 0.9314) methods. This indicates that the present invention can provide higher reconstruction accuracy and detail retention ability when the data is sufficient.

[0041] Under sparse sampling conditions (20 projections): In an extremely sparse scenario, the PSNRs of HWDSC and AWL1R are 26.59 dB and 26.38 dB respectively, and the SSIMs are 0.9191 and 0.9166 respectively, showing a significant improvement compared to FBP (PSNR = 14.50 dB, SSIM = 0.2805) and ART (PSNR = 22.20 dB, SSIM = 0.7395). This indicates that the present invention can still maintain the structural integrity and visual clarity of the image while significantly reducing the number of projections, reducing streak artifacts and noise interference.

[0042] Under intermediate sampling conditions: At 60 projections, the PSNRs of HWDSC and AWL1R are 29.30 dB and 29.69 dB respectively, and the SSIMs are 0.9619 and 0.9664 respectively, approaching the fully sampled effect, while the indicators of FBP and ART decrease significantly (22.72 dB, 0.5297 and 23.58 dB, 0.8907 respectively). This further verifies the superiority of the present invention under different sampling densities.

[0043] (2) Enhancing the robustness of the algorithm The present invention ensures the stability and adaptability of the algorithm under different numbers of projections and noise levels through multi-domain sparse representation (HWDSC) and adaptive regularization (AWL1R) strategies.

[0044] Adaptability to the change in the number of projections: Table 1 shows that when the number of projections is reduced from 180 to 20, the PSNR and SSIM of HWDSC and AWL1R decrease slightly. For example, at 30 projections, the SSIM of HWDSC and AWL1R are 0.9446 and 0.9443 respectively, while those of FBP and ART are only 0.3513 and 0.7815, indicating that the performance stability of the present invention under sparse conditions far exceeds that of traditional methods.

[0045] Noise suppression ability: Noise interference (such as 5% - 10% Gaussian noise) is added to the simulation data. HWDSC balances local details and global smoothness through multi - domain sparse constraints, and AWL1R dynamically adjusts the regularization strength through an adaptive threshold. The MSE of both under extremely sparse conditions (0.0022 and 0.0023 respectively at 20 projections) is much lower than that of FBP (0.0355) and ART (0.0060), showing a strong anti - noise ability.

[0046] (3)Support for low - dose imaging The present invention significantly reduces the radiation dose of CT scans by reducing the number of projections, while ensuring image quality and meeting the clinical requirements of low - dose imaging.

[0047] Reduction of radiation dose: Taking 60 projections as an example (a 66.7% reduction compared to 180), the PSNR of HWDSC and AWL1R are 29.30 dB and 29.69 dB respectively, and the SSIM are 0.9619 and 0.9664 respectively, approaching the reconstruction quality under full sampling (PSNR≈31 dB, SSIM≈0.975). In contrast, the PSNR of FBP and ART are only 22.72 dB and 23.58 dB respectively, and the SSIM drops significantly to 0.5297 and 0.8907. This indicates that the present invention can still provide high - quality images when the radiation dose is reduced by two - thirds, and is suitable for low - dose scans of sensitive populations such as children and pregnant women.

[0048] Extremely sparse scenarios: At 20 projections (an 88.9% reduction), the PSNR of HWDSC and AWL1R remains above 26 dB, and the SSIM is close to 0.92, far exceeding traditional methods, providing technical support for ultra - low - dose CT imaging.

[0049] (4)Optimization of computational efficiency The present invention uses the Fast Iterative Shrinkage Threshold Algorithm (FISTA) to optimize the reconstruction process, shortening the computational time while ensuring image quality. This benefits from the efficient iterative strategy of the FISTA algorithm and the reasonable step - size selection (based on Lipschitz constant estimation), meeting the real - time requirements of clinical and industrial applications.

[0050] (5)Wide application potential The present invention is not only applicable to medical CT imaging, but also can be extended to fields such as industrial non-destructive testing and micro-CT.

[0051] It demonstrates cross-field applicability: The data in Table 1 show that regardless of the number of projections, both HWDSC and AWL1R can provide stable and high-quality reconstruction results. For example, at 90 projections, the PSNR reaches 29.65 dB and 30.20 dB respectively, and the SSIM is 0.9647 and 0.9714, which can be applicable to the detection of complex structures in industrial CT and high-resolution imaging in biological science research.

[0052] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, and the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the specified functions in Figure 1 one process or multiple processes and / or blocks Figure 1 one block or multiple blocks.

[0053] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including instruction means, and the instruction means implements the specified functions in Figure 1 one process or multiple processes and / or blocks Figure 1 one block or multiple blocks.

[0054] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable device provide steps for implementing the specified functions in Figure 1 one process or multiple processes and / or blocks Figure 1 one block or multiple blocks.

[0055] Specific embodiments are applied in the present invention to elaborate on the principles and implementation manners of the present invention. The description of the above embodiments is only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, based on the idea of the present invention, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation to the present invention.

[0056] Those of ordinary skill in the art will realize that the embodiments described herein are provided to assist the reader in understanding the principles of the present invention, and it should be understood that the scope of protection of the present invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific deformations and combinations that do not depart from the essence of the present invention based on these technical revelations disclosed in the present invention, and these deformations and combinations are still within the scope of protection of the present invention.

Claims

1. A sparse angle CT image reconstruction method based on the theory of compressive sensing, characterized in that, It includes the following steps: S1. Load the projection data and initialize the parameter settings, where the parameters include the number of projections, the number of iterations, and the tolerance; S2. Apply the FBP method to the projection data and output the initial reconstructed image; S3. Input the obtained initial reconstructed image into the ART reconstruction module and perform ART module iteration; S4. Perform external optimization on the image output after ART module iteration, and update the projection image after the optimization ends.

2. The sparse-angle CT image reconstruction method based on the compressive sensing theory according to claim 1, wherein The specific steps of S3 are as follows: S31. Input the obtained initial reconstructed image into the ART reconstruction module and save the current image; S32. Traverse all the number of angles, extract the current angle and calculate the forward projection of the current angle, calculate and save the projection error of the current angle, and end the inner loop of the angle; S33. Perform back projection on the projection error at the current angle and update the current image, perform FISTA denoising and update the current image; S34. Determine whether the updated image converges. If it does not converge, continue to perform ART module iteration until it converges.

3. The sparse-angle CT image reconstruction method based on the compressive sensing theory according to claim 2, wherein The specific formula of the ART module in S34 is: where x is the image to be reconstructed, a i is the i-th row of the projection matrix A, y i is the projection data, λ is the relaxation parameter, k is the number of iterations, T is the matrix transpose, is the norm calculation.

4. A sparse angle CT image reconstruction method based on the theory of compressive sensing according to claim 1, characterized in that In S4, the external optimization adopts an optimization algorithm based on hybrid wavelet and DCT coefficient constraint or an optimization algorithm based on adaptive wavelet transform and L1 regularization.

5. A sparse angle CT image reconstruction method based on the theory of compressive sensing according to claim 4, characterized in that, The optimization algorithm based on hybrid wavelet and DCT coefficient constraint uses the reconstructed image of filtered back projection as the initial value, and its optimization objective function is expressed as: Where, Ψ x is the wavelet transform operator, λ1 and λ2 are constants, x is the image to be reconstructed, DCT is the discrete cosine transform operator, is the norm calculation, y is the measurement data, and A is the projection matrix.

6. A sparse angle CT image reconstruction method based on the theory of compressive sensing according to claim 4, characterized in that The optimization algorithm based on adaptive wavelet transform and L1 regularization uses the reconstructed image of filtered back projection as the initial value, and its optimization objective function is expressed as: Where x is the image to be reconstructed, μ is the L1 regularization parameter in the image domain, is the norm calculation, y is the measurement data, and A is the projection matrix; λ is the regularization parameter in the wavelet domain, and its adaptive threshold adjustment method is as follows: ； In the formula, a is a constant, n is the current number of projections, N is the maximum number of projections, and i and I are the current and total number of iterations respectively.

Images