CT reconstruction system and method based on discrete Gaussian representation

Abstract

The invention belongs to the technical field of medical image reconstruction, and particularly relates to a CT reconstruction system and method based on discrete Gaussian representation. According to the technical scheme, the method comprises the steps that a Gaussian function set is constructed, and a three-dimensional voxel space is directly represented through learnable Gaussian parameters; carrying out parallel aggregation on Gaussian function contribution values to a three-dimensional voxel grid by adopting a discretization volume reconstruction technology; establishing a geometric projection transformation model to reversely project the three-dimensional reconstruction volume to a two-dimensional projection domain, and minimizing the difference of the projection domain through joint optimization of Gaussian parameters; and introducing an adaptive density control mechanism to dynamically adjust Gaussian distribution. Compared with an existing method, the method has the advantages that an end-to-end discretization Gaussian function representation mode is adopted, global optimization of any CT geometric configuration is supported, high-quality rapid reconstruction within 5-15 minutes is achieved while the detail integrity of the anatomical structure is kept, and a high-precision and high-robustness three-dimensional reconstruction solution is provided for low-dose CT scanning and emergency treatment scenes.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and particularly relates to a three-dimensional voxel reconstruction system and method for computed tomography (CT) based on discretized Gaussian representation, which is particularly suitable for fast and high-precision three-dimensional volume reconstruction under sparse projection data. Background Art

[0002] The three-dimensional reconstruction technology of computed tomography (CT) is the cornerstone of clinical image diagnosis, and its core lies in inverting the voxel spatial distribution from multi-angle X-ray projection data. The current mainstream methods can be divided into two categories: reconstruction technologies based on deep learning and reconstruction technologies based on single-instance optimization. However, both of them have essential technical defects and are difficult to meet the comprehensive clinical requirements for reconstruction speed, cross-device generalization ability, and anatomical fidelity.

[0003] Although the traditional Filtered Back Projection (FBP) algorithm has the ability of real-time reconstruction, it will produce serious noise and streak artifacts under low-dose conditions; the Iterative Reconstruction (IR) algorithm can suppress noise through mathematical model optimization, but it has the inherent problems of taking up to dozens of minutes in calculation time and tissue texture distortion caused by over-smoothing. In recent years, the emerging Deep Learning Reconstruction (DLR) uses convolutional neural networks to learn the projection-voxel mapping relationship from a large amount of labeled data. Although it has made progress in noise suppression, its supervised training mechanism relying on tens of thousands of pairs of paired data (such as low-dose and conventional-dose images) makes the model severely limited by specific scanning parameters (such as fan-beam CT geometry) or anatomical parts (such as chest scans). When migrated to cone-beam CT or head scans, it will lead to a significant performance decline, forcing medical institutions to retrain the model for different scenarios, significantly pushing up the deployment cost.

[0004] To avoid the defect of data dependence, the method based on instance optimization realizes reconstruction through iterative optimization of a single sample. Typical representatives include Neural Radiance Fields (NeRF) and 3D Gaussian Splatting (3DGS) technology. However, NeRF-like methods require hours of training time to fit the radiation field parameters and cannot meet the minute-level response requirements of emergency scenarios; although the 3DGS method improves the rendering speed through discretized Gaussian basis element modeling (such as the adaptive rasterization strategy adopted in Chinese Patent CN119600138A), the essential conflict between its continuous representation and the CT discrete voxel grid leads to reconstruction accuracy errors, and the rasterization process increases the time and video memory overhead.

[0005] The above dilemmas reveal three core contradictions in current technologies: First, the training cost of deep learning reconstruction models is too high and the generalization ability is insufficient; Second, the method based on instance optimization is not compatible with the mathematical representation of discrete voxel grids, resulting in performance loss, and the non-direct connection between rasterization rendering and three-dimensional volume reconstruction limits the reconstruction efficiency; Third, there is an urgent clinical need for end-to-end three-dimensional reconstruction to be completed within minutes, but there is always a trade-off between reconstruction efficiency and accuracy. These bottlenecks severely restrict the deepening of the application of CT technology in real-time sensitive scenarios such as emergency medicine, and there is an urgent need for a CT reconstruction method that can achieve high precision, high efficiency, and high generalization at the same time. Summary of the Invention

[0006] In view of the technical bottlenecks in computer tomography (CT) image reconstruction, the present invention proposes a CT reconstruction system and method based on discrete Gaussian representation (Discretized Gaussian Reconstruction, abbreviated as DGR). Through a modular system architecture, combined with innovative Gaussian function representation, initialization strategy, reconstruction algorithm, and adaptive density optimization method, efficient, flexible, and high-precision image reconstruction is achieved. The specific technical solutions are as follows:

[0007] 1. System Architecture

[0008] The present invention designs a modular system architecture that supports the complete processing flow from CT projection data to high-quality three-dimensional reconstructed images. The system includes the following modules:

[0009] Input Module: Used to receive CT projection data with different geometric structures, including but not limited to cone beam (Cone Beam) scanning mode, fan beam (Fan Beam) scanning mode, etc.;

[0010] Initialization Module: Generate the initial parameters of the Gaussian function based on the results of the FDK (Feldkamp-Davis-Kress) algorithm;

[0011] Discrete Reconstruction Module: Combine multiple Gaussian functions into a three-dimensional volume through a discretized volume reconstruction method;

[0012] Parameter Update Module: Optimize the Gaussian function parameters in the projection domain to minimize the error between the predicted projection and the input projection;

[0013] Density Optimization Module: Dynamically adjust the density and distribution of the Gaussian function according to the reconstruction quality to balance accuracy and efficiency;

[0014] Output Module: Present the reconstructed three-dimensional volume in a visual manner.

[0015] Among them, the discrete reconstruction module, the parameter update module, and the density optimization module constitute a multi-round iterative update unit, and the reconstruction performance is improved by cyclically executing parameter optimization. The system execution process includes:

[0016] a) Initialization stage: Sequentially activate the input module and the initialization module to establish the initial parameters of the Gaussian function;

[0017] b) Multi-round iterative optimization stage: Cyclically execute the basic iteration including the discrete reconstruction module and the parameter update module, and activate the density optimization module every preset round until the convergence condition is met;

[0018] c) Result output stage: Generate a visual three-dimensional reconstruction result through the output module.

[0019] 2. Three-dimensional representation of the Gaussian function

[0020] To solve the problem of sparse data in low-dose CT images, the present invention uses a Gaussian function to parametrically represent the three-dimensional volume. The volume data V is represented as a weighted sum of a set of Gaussian functions:

[0021]

[0022] Among them:

[0023]

[0024] μ i is the center point of the Gaussian function, Σ i is the covariance matrix, and I i is the intensity of the Gaussian function.

[0025] This representation has the following advantages:

[0026] It can efficiently describe the sparsity of low-dose CT data;

[0027] It is convenient for subsequent fast reconstruction and optimization operations;

[0028] It can adapt to multiple projection perspectives and angles without re-designing the network architecture.

[0029] 3. Efficient initialization strategy

[0030] To improve the reconstruction efficiency, the present invention designs a Gaussian function initialization strategy based on the FDK result.

[0031] First, calculate the two-norm of the gradient of the volume data V FBP generated by filtered back projection:

[0032]

[0033] Then, calculate V gnormThe median value of non - negative values, denoted as M:

[0034]

[0035] Next, select V FBP k points in V whose gradient values are closest to M as the centers of the k initialized Gaussian functions;

[0036] The covariance matrix Σ of each Gaussian function i is initialized through its geometric relationship with adjacent Gaussian distributions. The steps are as follows:

[0037]

[0038] where N k (i) represents the set of indices of the k Gaussian functions closest to μ i nearest.

[0039] Finally, initialize the covariance matrix as a diagonal matrix:

[0040] Σ i = diag(σ i , σ i , σ i )·γ (6)

[0041] The value range of the scaling factor γ is recommended to be [0.1, 0.5], which is used to adjust the initial coverage range of the Gaussian distribution.

[0042] 4. Three - dimensional volume reconstruction method

[0043] The present invention adopts a direct three - dimensional reconstruction method to generate a three - dimensional volume through Gaussian functions. The specific process is as follows:

[0044] Divide the three - dimensional space into discrete voxel grids, and the intensity at each voxel position is determined by the cumulative contribution of all Gaussian functions; the intensity calculation formula for each voxel V p is:

[0045]

[0046] where x p is the central coordinate of the voxel, and G(x p , μ i , Σ i ) is the contribution of the i - th Gaussian function to x p ;

[0047] To improve the calculation efficiency, only consider the contribution of Gaussian functions within the local area to the voxel, and ignore the Gaussian functions far from the voxel center. Through this direct three - dimensional reconstruction method, the present invention can generate a high - fidelity three - dimensional volume under low - dose data conditions.

[0048] 5. Discretized Volume Reconstruction Strategy

[0049] To further accelerate the reconstruction process, the present invention designs a local restriction and matrix decomposition strategy. The specific operations include:

[0050] Restrict the action range of the Gaussian function within its neighboring grid area, and define the local influence range as a cuboid. The number of voxels of the cuboid in the width, height, and depth directions are denoted as w0, h0, and c0 respectively.

[0051] To improve the reconstruction efficiency and reduce the time and space complexity, the present invention proposes a matrix decomposition method applicable to fast volume reconstruction, which decomposes high-dimensional operations into multiple low-dimensional operations to significantly reduce the computational complexity.

[0052] Denote the difference between the center μ of the continuous Gaussian function and the center of the discrete Gaussian function as Δμ, that is:

[0053]

[0054] where represents the discretized rounding operation on the continuous center coordinate μ.

[0055] To intuitively describe the transformation of tensor dimensions, the following formula follows the Einstein summation convention, where the subscripts represent tensor dimensions. First, calculate the quadratic form of the inverse of the covariance matrix and the local grid deviation:

[0056]

[0057] Then calculate the linear term of the deviation and the center coordinate offset:

[0058]

[0059] Next, calculate the constant term of the center coordinate offset:

[0060]

[0061] where B is the deviation between the local grid and the Gaussian center, Σ -1 is the inverse of the covariance matrix, and w0, h0, c0 respectively represent the number of voxels of the truncated region in the width, height, and depth directions.

[0062] Then calculate the Mahalanobis distance of the discretized Gaussian distribution:

[0063]

[0064] where the subscript n represents the nth Gaussian function, and w0, h0, c0 are the local cuboid sizes;

[0065] The contribution of all Gaussian functions within the moment body can be characterized as follows:

[0066]

[0067] where I is the tensor composed of the intensities of all Gaussian functions in formula (1);

[0068] Finally, the contributions Γ of all the above Gaussian functions are superimposed into the three-dimensional voxel grid V, i.e., the discretized volume reconstruction is achieved.

[0069] 6. Projection Domain Parameter Update Strategy

[0070] The present invention defines a parameter update method in the projection domain to achieve iterative optimization of Gaussian functions.

[0071] First, the three-dimensional volume is projected into the projection domain through the Radon transform:

[0072] S pred = Radon(V pred ) (15)

[0073] where the Radon transform here adopts the same Radon transform method as the original CT projection data, and the generated three-dimensional volume is converted into the projection domain, so that for any CT sampling mode, such as cone-beam CT and fan-beam CT, three-dimensional volume reconstruction can be completed without changing the method structure.

[0074] Subsequently, the error between the predicted projection (S pred ) and the input projection (S input ) is minimized. Denote the loss function as L, and we have:

[0075] Loss = L(S pred , S input ) (16)

[0076] where the core idea of this method is to calculate the loss in the projection domain to achieve iterative optimization, so that data in other three-dimensional volume domains do not need to be used, thereby reducing the radiation dose borne by patients and reducing the computational overhead. Therefore, there are many different definition methods for L, as long as the error between calculating S pred and S input is satisfied. For example, the L2 loss is used:

[0077]

[0078] To further improve the performance, the Total Variation loss and the SSIM loss can be added. Their definitions are as follows:

[0079]

[0080] Among them, the TV loss is directly calculated for the reconstructed three-dimensional volume V pred by substituting each voxel into the above formula. It should be noted that the supervision signal of the three-dimensional volume is not used here, so calculating the loss does not conflict with the projection domain optimization.

[0081]

[0082] Among them, the SSIM loss is calculated using S pred and S input by substituting into the above formula.

[0083] Then the total loss is:

[0084] L total = λ1L L2 + λ2L SSIM + λ3L TV (20)

[0085] Among them, λ1, λ2, λ3 are weight hyperparameters used to balance the contributions of different loss terms;

[0086] The above is only an example of calculating one type of loss. For the calculation of other types of loss functions, as long as the optimization in the projection domain is satisfied, the projection domain parameter update strategy of the present invention can also be realized.

[0087] When the number of iteration rounds reaches the set maximum round limit, or L total is lower than the loss threshold ∈ i , the iterative optimization step terminates.

[0088] 7. Adaptive density optimization strategy

[0089] The present invention uses an adaptive density optimization strategy based on 3DGS to dynamically adjust the distribution and density of Gaussian functions through cloning, splitting, and pruning operations; every α rounds (where α is a hyperparameter), the density optimization module is executed once;

[0090] The density optimization module can dynamically adjust the distribution and density of Gaussian functions, and its specific operations include:

[0091] (1) Cloning: For the Gaussian functions in the low-intensity region, clone them into η1 copies. The central coordinates of each cloned Gaussian function remain unchanged, and the intensity decays to the original

[0092] μ clone = μ original (21)

[0093]

[0094] Among them, μoriginal is the central coordinate of the Gaussian function to be cloned, μ clone is the central coordinate of the cloned Gaussian function; I original is the intensity of the Gaussian function to be cloned, I clone is the intensity of the cloned Gaussian function; the cloned Gaussian function shares the covariance matrix with the original function;

[0095] (2) Splitting: Split the Gaussian function in the reconstructed complex region into η2 parts. The central coordinates of each split Gaussian function remain unchanged, and the covariance decays to the original

[0096] μ split = μ original (23)

[0097]

[0098] where, μ original is the central coordinate of the Gaussian function to be split, μ split is the central coordinate of the split Gaussian function, Σ original is the covariance of the Gaussian function to be split, Σ split is the covariance of the split Gaussian function;

[0099] (3) Pruning: Remove redundant Gaussian functions that contribute less to the loss or have gradients close to zero. Specifically, remove I i below the threshold ∈ I or the gradient norm continuously below the threshold ∈ g of the Gaussian function during the parameter update process, where θ i represents μ i , Σ i or I i . This strategy can dynamically balance the computational efficiency and the reconstruction accuracy to meet the requirements of different clinical scenarios.

[0100] Compared with the prior art, the advantages of the present invention are reflected in:

[0101] By constructing a complete CT reconstruction method system from underlying representation modeling to reconstruction algorithm implementation and then to optimization strategy design, the technical bottlenecks existing in traditional methods are effectively overcome. The discretized Gaussian basis function representation technology is adopted to make the reconstruction objective function form a mathematical isomorphism with the discretized nature of the three-dimensional volume, and it is significantly superior to existing mainstream CT reconstruction methods in terms of image quality indicators. Specifically: compared with deep learning reconstruction methods, it still maintains a 1dB PSNR advantage under unsupervised conditions on the AAPM-Mayo LDCT public dataset; compared with the most advanced instance optimization reconstruction methods, it achieves a 1-2dB PSNR performance improvement on the FIPS dataset. In terms of computational efficiency, the discretized volume reconstruction framework proposed in the present invention can complete high-precision three-dimensional reconstruction within 5-15 minutes through algorithm architecture innovation, achieving a 3-5 times improvement in computational efficiency compared with existing methods. The innovative advantages of the present invention can be further summarized into the following four dimensions:

[0102] Accuracy: The quality of the reconstructed image is significantly improved through the adaptive density optimization algorithm for the discretized representation space;

[0103] Efficiency: The computational complexity is greatly reduced through discretized volume reconstruction and local restriction strategies;

[0104] Adaptability: It does not rely on the supervised training of full-dose CT data and supports diverse projection angle configurations;

[0105] Generalization: Based on the directional reconstruction mechanism for patient-specific data, it can adaptively match the imaging requirements of different anatomical structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0106] Figure 1 It is a schematic diagram of the system architecture of the present invention

[0107] Figure 2 It is a three-dimensional representation diagram of the Gaussian function

[0108] Figure 3 It is a schematic representation diagram of the discretized Gaussian function

[0109] Figure 4 It is a schematic diagram of the discretized volume reconstruction algorithm of the present invention

[0110] Figure 5 It is a schematic diagram of the result of iterative reconstruction optimization of the present invention

[0111] Figure 6 It is a cross-sectional comparison diagram of the experimental results of the embodiment of the present invention

[0112] Figure 7 It is a three-dimensional comparison diagram of the experimental results of the embodiment of the present invention DETAILED DESCRIPTION OF THE INVENTION

[0113] To clearly illustrate the technical solution of the present invention, the following describes the specific implementation manner of the present invention in detail in combination with the implementation steps. It should be noted that the following embodiments are used to further illustrate the technical solution of the present invention, but the protection scope of the present invention is not limited to these embodiments. Equivalent replacements, modifications or reasonable expansions made by those skilled in the art within the scope defined by the claims, based on the core idea of the present invention, all belong to the protection scope of the present invention.

[0114] 1. System architecture design

[0115] The CT system of this embodiment includes an input module, an initialization module, a discrete reconstruction module, a parameter update module, a density optimization module, and an output module.

[0116] 2. System input

[0117] (1) The input module receives the projection data S generated by CT scanning input ;

[0118] (2) Convert the input projection data into a format adapted to the Gaussian representation to support subsequent calculations.

[0119] 3. Gaussian function initialization

[0120] The initialization module extracts the initial parameters of the Gaussian function according to the volume data V generated by FDK FDK , including the Gaussian center point μ, the covariance matrix Σ, and the intensity I. The specific steps are as follows:

[0121] (1) Calculate the two-norm of the gradient of V FDK :

[0122]

[0123] (2) Calculate the median M according to the gradient norm:

[0124]

[0125] Select the k points in V FBP whose gradient values are closest to M as the initialization Gaussian function center points. In this embodiment, k = 150000 is taken; (3) Initialize the covariance matrix based on the average value of the squared Euclidean distances between adjacent Gaussian points:

[0126]

[0127] Σ i = diag(σ i , σ i , σ i )·γ (28)

[0128] where N k(i) represents the k nearest neighbor points of the i-th Gaussian point. In this embodiment, γ = 0.1 is taken;

[0129] (4) The initial value of the Gaussian intensity I is set to the intensity of the corresponding voxel in the FDK reconstructed volume:

[0130] I i = V FDK (μ i ) (29) 4. Use discretized volume reconstruction to complete the fast reconstruction of the Gaussian function into a three-dimensional volume. The specific steps are as follows:

[0131] (1) Limit the action range of the Gaussian function to a local area and define a rectangular box B:

[0132] B = {(x, y, z) ∣ |x - μ x | ≤ w0, |y - μ y | ≤ h0, |z - μ z | ≤ c0} (30)

[0133] where w0, h0, c0 are the sizes of the local range; in this embodiment, w0 = 17, h0 = 17, c0 = 17;

[0134] (2) According to formulas (9) - (13), use matrix decomposition to optimize the calculation of local contributions:

[0135] D 2 = B T Σ -1 B (31)

[0136] (3) Superimpose the calculated contributions of the Gaussian function onto the three-dimensional voxel grid to obtain the reconstructed volume.

[0137] 5. Optimize the parameters μ, Σ, and I of the Gaussian function by minimizing the loss function in the projection domain. The specific steps are as follows: (1) Project the three-dimensional volume onto the projection domain through the Radon transform:

[0138] S pred = Radon(V pred ) (32)

[0139] where V pred is the three-dimensional volume represented by the current Gaussian function;

[0140] (2) Define the projection domain loss function. In this embodiment, the simplest L2 loss is adopted:

[0141]

[0142] where S input is the input projection, S predis a predictive projection;

[0143] In this embodiment, the upper limit of the maximum number of rounds is set to 1000, and the loss threshold ∈ i = 0.001.

[0144] (3) Use the gradient descent method to optimize the loss function and update the Gaussian parameters.

[0145] 6. Adaptive density optimization

[0146] The density optimization module dynamically adjusts the distribution and density of the Gaussian function. According to formulas (21)-(24), the following density optimization operations are performed on the Gaussian function. In this embodiment, α = 100, η1 = 2, η2 = 2, ∈ g = 1e-5:

[0147] (1) Cloning: For regions with insufficient reconstruction accuracy, the original Gaussian function is cloned into two Gaussian functions, and at the same time, its intensity is reduced to keep the total intensity constant;

[0148] (2) Splitting: In regions with complex local reconstructions, a single Gaussian function is split into two sub-Gaussian functions with smaller covariance to more finely describe the structural details;

[0149] (3) Pruning: Redundant Gaussian functions with gradients close to zero or small contributions are pruned, thereby reducing the computational burden and avoiding overfitting.

[0150] 7. Output and display

[0151] The finally reconstructed three-dimensional volume V is displayed through the output module, supporting multi-view switching and interactive analysis, including cross-sectional view display (as shown in the appendix Figure 6 ), and three-dimensional volume display (as shown in the appendix Figure 7 ).

[0152] 8. Table 1 shows the experimental results of this implementation case on the AAPM-Mayo LDCT dataset. This verification uses a sparse fan-beam CT scanning mode, and 10 abdominal CT clinical cases are selected as test examples. The experimental comparison scheme includes advanced DLR (Deep Learning Reconstruction) methods based on deep learning and other representative algorithms. The literature sources of each comparison method are marked behind the corresponding method names. To fully verify the effectiveness of the method, the following control dimensions are specifically set: 1) In terms of the scale of training data, the capacity of the training set for DLR-based methods is clearly marked in the "Training Data" column; 2) In the dimension of the number of projections, four sparse sampling schemes of "180-view" (representing reconstruction using 180 projection views), "120-view", "90-view", and "60-view" are set. 3) At the level of performance evaluation, the experiment uses two core indicators, peak signal-to-noise ratio (PSNR / dB) and structural similarity index (SSIM), for quantitative evaluation. The optimal performance data under each indicator dimension are highlighted in bold font. It should be noted that all comparison algorithms are strictly executed according to the parameter configurations and implementation conditions recorded in their original literature.

[0153] Table 1 Comparison of Reconstruction Performance with DLR-based Methods on the AAPM-Mayo LDCT Dataset (PSNR / SSIM)

[0154]

[0155] In the reconstruction performance verification of the embodiment of the present invention on the AAPM-Mayo LDCT standard dataset, a technical breakthrough of zero additional training data requirement is innovatively achieved. As shown in Table 1, this method is significantly superior to the CT reconstruction methods disclosed in internationally authoritative academic publications (such as journals / conferences like IEEE Transactions on Medical Imaging) cited by the comparison group in terms of reconstruction accuracy (quantitative indicators include PSNR and SSIM). This comparative experiment strictly follows the technical parameter settings of the original literature, fully demonstrating the advantages of high precision and high generalization of the present invention.

[0156] 9. Table 2 shows the experimental results of this embodiment on the FIPS dataset. This verification experiment is carried out based on the FIPS standard dataset, adopting a sparse cone-beam CT scanning protocol, and selecting 15 CT scanning instances with different object morphologies as the test benchmarks. The experimental comparison scheme includes an advanced algorithm system optimized based on instances, and the technical sources of each comparison method are marked in parentheses behind the corresponding method names. This embodiment specifically sets two comparison experiments with different training intensities: DGR(iter = 300) corresponds to the conventional training mode with 300 iterations, and DGR(iter = 1000) corresponds to the intensive training mode with 1000 iterations. The experimental settings include the following core verification dimensions: 1) At the level of the number of projections, three levels of sparse sampling are set, namely "75-view" (representing the reconstruction benchmark of 75 projection views), "50-view", and "25-view"; 2) At the level of performance evaluation, the experiment uses three core indicators, namely peak signal-to-noise ratio (PSNR / dB), structural similarity index (SSIM), and algorithm running time, for quantitative evaluation. The optimal performance data under each indicator dimension are highlighted in bold font. It should be specifically stated that all comparison algorithms strictly reproduce the technical solutions of their original literatures.

[0157] Table 2 Performance Comparison with Instance-Optimized Reconstruction Methods under the FIPS Dataset (PSNR / SSIM / Time)

[0158]

[0159] As shown in the experimental data in Table 2, in the comparison with the advanced method system optimized based on instances, the embodiment of the present invention shows double technical breakthroughs in accuracy and efficiency, and can complete high-precision CT reconstruction within a time window of 5 - 15 minutes (the specific time consumption is positively correlated with the number of scanning views), further corroborating the technical advantages of high precision, high efficiency, and high generalization of the present invention, and supporting the minute-level response in the emergency scenario.

[0160] The specific embodiments of the present invention have been described in detail above. However, it should be emphasized that these embodiments are only for illustrative purposes and not for limiting the scope of protection. The scope of protection of the present invention shall be subject to the claims, and the description and drawings can be used to explain the claims. Equivalent technical solutions obtained by those skilled in the art through conventional experiments, logical reasoning, or limited substitutions on the basis of understanding the technical solutions of the present invention and in combination with the prior art, as long as they do not deviate from the core features defined by the claims of the present invention, all fall within the scope of protection of the present invention.

[0161] References

[0162] [1]Kyong Hwan Jin,Michael T McCann,Emmanuel Froustey,and Michael Unser.Deep convolutional neural network for inverse problems in imaging.IEEE transactions on image processing,26(9):4509–4522,2017.

[0163] [2]Hoyeon Lee,Jongha Lee,Hyeongseok Kim,Byungchul Cho,and Seungryong Cho.Deep-neural-network-based sinogram synthesis for sparse-view ct image reconstruction.IEEE Transactions on Radiation and Plasma Medical Sciences,3(2):109–119,2018.

[0164] [3]Liutao Yang,Rongjun Ge,Shichang Feng,and Daoqiang Zhang.Learning projection views for sparse-view ct reconstruction.In Proceedings of the 30th ACM International Conference on Multimedia,pages 2645–2653,2022.

[0165] [4]Bing Guan,Cailian Yang,Liu Zhang,Shanzhou Niu,Minghui Zhang,Yuhao Wang,Weiwen Wu,and Qiegen Liu.Generative modeling in sinogram domain for sparse-view ct reconstruction.IEEE Transactions on Radiation and Plasma Medical Sciences,2023.

[0166] [5] Kai Xu, Shiyu Lu, Bin Huang, Weiwen Wu, and Qiegen Liu. Stage-by-stage wavelet optimization refinement diffusion model for sparse-view ctreconstruction. IEEE Transactions on Medical Imaging, 2024.

[0167] [6] Guangming Zang, Ramzi Idoughi, Rui Li, Peter Wonka, and Wolfgang Heidrich. Intratomo: self-supervised learning based tomography via sinogram synthesis and prediction. In Proceedings of the IEEE / CVF International Conference on Computer Vision, pages 1960–1970, 2021.

[0168] [7] Ruyi Zha, Yanhao Zhang, and Hongdong Li. Naf: neural attenuation fields for sparse-view cbct reconstruction. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 442–452. Springer, 2022.

[0169] [8] Yuanhao Cai, Jiahao Wang, Alan Yuille, Zongwei Zhou, and Angtian Wang. Structure-aware sparse-view x-ray 3d reconstruction. In Proceedings of the IEEE / CVF Conference on Computer Vision and Pattern Recognition, pages 11174–11183, 2024.

[0170] [9]Ruyi Zha, Tao Jun Lin, Yuanhao Cai, Jiwen Cao, Yanhao Zhang, and Hongdong Li. R2-gaussian: Rectifying radiative gaussian splatting for tomographic reconstruction. In Advances in Neural Information Processing Systems (NeurIPS), 2024.

Claims

1. A three-dimensional volume reconstruction system for computed tomography (CT) based on discrete Gaussian representation, characterized in that, The system includes the following modules: An input module for receiving CT projection data with different geometric structures; An initialization module for generating an initial parameter set of Gaussian functions based on the filtered back-projection (FDK) algorithm; A discrete reconstruction module for combining multiple Gaussian functions through a discretized volume reconstruction method to form a three-dimensional volume model; A parameter update module for optimizing Gaussian function parameters in the projection domain to minimize the error between the predicted projection and the reference projection; A density optimization module for dynamically adjusting the spatial density distribution of Gaussian functions according to the reconstruction quality; An output module for converting the reconstructed three-dimensional volume data into a visual expression form.

2. The system according to claim 1, wherein The initialization module performs the following operations: Calculate the gradient two-norm of the volume data V generated by filtered backprojection FDK ​ Calculation Median of non-zero values Select V FDK Select the k points with the nearest M in the middle gradient two-norm as the initial center point μ of the Gaussian function i ; Initialize the covariance matrix Σ according to the Euclidean distance between adjacent Gaussian center points i , specifically And let Σ i = diag(σ i , σ i , σ i )·γ, where γ is an adjustable scaling factor.

3. The system according to claim 1, wherein The discrete reconstruction module performs the following operations: Dividing the three-dimensional space into a discrete voxel grid and defining the scope of each Gaussian function as a rectangular parallelepiped region of w0×h0×c0; Calculate the discretized Mahalanobis distance within the rectangular region And calculate the total contribution of all Gaussian functions: Performing superposition integration of the total contribution Γ of the above Gaussian functions in the three-dimensional voxel grid to complete the discretized volume reconstruction.

4. The system according to claim 1, wherein The parameter update module: Configuring a multi-geometric compatibility mechanism that can adaptively process the input of cone-beam CT and fan-beam CT data with different numbers of projections without reconstructing the system architecture; Constructed based on the Radon transform parameters corresponding to the original CT projection data, converting the three-dimensional volume data Vpred to the projection domain to generate the predicted projection data Spred = Radon(Vpred); Construct a projection domain error calculation model to implement the loss function calculation between the predicted projection data Spred and the input projection data S input ; establish a gradient backpropagation channel, and iteratively update the center coordinate μ of the Gaussian function through the backpropagation algorithm i , covariance matrix Σ i , and intensity coefficient I i .

5. The system according to claim 1, characterized in that, The density optimization module performs the following operations: Cloning operation: Copy the Gaussian function in the low-density region into η1 copies, keeping the center position and covariance matrix unchanged, and adjusting the intensity to Splitting operation: Split the Gaussian function in the high-density region into η2 parts, keeping the center position and the total intensity unchanged, and reducing the covariance matrix to Pruning operation: Remove Gaussian functions where the intensity is below the threshold ∈ I or the gradient norm continuously below the threshold ∈ g , where θ i represents μ i , Σ i or I i parameters.

6. The system according to claim 1, characterized in that The execution logic of each module is as follows: Executing the input module and the initialization module once to obtain the initial Gaussian parameter set from the projection data; Perform multiple rounds of iterative updates. In each round, execute the discrete reconstruction module and the parameter update module, and execute the density optimization module every α rounds; when the number of iterations reaches the preset upper limit or the total loss L total is lower than the threshold ∈ i terminate the optimization Finally, presenting the reconstruction result in the form of a two-dimensional cross-section and three-dimensional volume rendering through the output module.

7. A CT three-dimensional volume reconstruction method based on discrete Gaussian representation, characterized in that, Including: Performing Gaussian parameterization characterization on the three-dimensional volume; Decomposing the Gaussian contribution calculation into four tensor sub-operations through a matrix decomposition strategy; Adopting a composite loss function of linear combination during projection domain optimization; Dynamically adjusting the spatial distribution density of Gaussian functions according to the reconstruction quality.

8. The method according to claim 7, wherein The parametric representation of the Gaussian function includes the center point μ of the Gaussian function i , the covariance matrix Σ i , the intensity I of the Gaussian function i , and the three-dimensional volume V is expressed as: where the Gaussian function G(x, μ i , Σ i ) is defined as:

9. The method according to claim 7, wherein The discretized volume reconstruction uses a matrix decomposition strategy to decompose the Gaussian function contribution calculation into four sub-tensor operations to save computational overhead, specifically including: Calculating the quadratic form of the inverse of the covariance matrix and the local grid deviation; Calculating the linear term of the deviation and the center coordinate offset; Calculating the constant term of the center coordinate offset; Combining to form the Mahalanobis distance of the discretized Gaussian distribution; 10. The method according to claim 7, characterized in that The construction method of the composite loss function adopted during projection domain optimization includes: Establish a multi-loss collaborative optimization mechanism to decompose the difference calculation between the prediction projection Spred and the input projection S input into multiple loss terms, including L2 loss, Total Variation loss, and SSIM loss; The composite loss function L is constructed by weighted combination of the loss terms with hyperparameters total = λ1L L2 + λ2L SSIM + λ3L TV , where λ1, λ2, and λ3 are adjustable hyperparameter weights

Images