Image reconstruction method and system based on sparse graph set

Abstract

An image reconstruction method and system based on a sparse graph set. The method comprises the following steps: using CT scanning to obtain a sparse angle projection graph set; performing back projection on the sparse angle projection graph set to obtain an initial iteration value of volume data; calculating initial values of three directional gradient fields of the initial iteration value of the volume data; on the basis of the initial iteration value of the volume data and the initial values of the three gradient fields, respectively calculating orthographic projection and residual back projection; performing, by means of a neural network and by using the initial values of the three gradient fields, regularization processing on data having undergone orthographic projection and data having undergone residual back projection; performing addition on an integral of the volume data and the gradient fields to serve as an initial value of next iteration, iterating again until completion and outputting a result. The method can obtain projection data required by reconstruction in a short period of time, and can effectively process the sparse angle / finite angle problem to obtain a high-quality image reconstruction result.

Description

Image reconstruction method and system based on sparse atlas Technical Field

[0001] The present invention relates to an image reconstruction method based on a sparse atlas, and also relates to a corresponding image reconstruction system, belonging to the technical field of digital image processing. Background Art

[0002] When a static CT system scans an object, its radiation source and detector remain stationary or rotate only within a small range. The technical advantages of this scanning method include: first, it no longer relies on a complex slip ring structure, simplifying the equipment design; second, due to the reduction of mechanical movement, the scanning speed is significantly improved, thereby improving the temporal resolution; in addition, the image smear effect that may be caused by high-speed movement is avoided, thereby improving the spatial resolution. However, since the projection geometry of the static CT system is different from that of the traditional spiral CT system, this requires the development of new iterative reconstruction algorithms to effectively process projection data with sparse or limited angles. At the same time, it is also necessary to develop efficient forward / backward projection operators suitable for multi-flat panel detector stitching to meet the imaging needs of static CT systems.

[0003] Summary of the Invention

[0004] The primary technical problem to be solved by the present invention is to provide an image reconstruction method based on a sparse atlas.

[0005] Another technical problem to be solved by the present invention is to provide an image reconstruction system based on a sparse atlas.

[0006] In order to achieve the above technical objectives, the present invention adopts the following technical solutions:

[0007] According to a first aspect of an embodiment of the present invention, a sparse atlas-based image reconstruction method is provided, comprising the following steps:

[0008] Using CT scanning, a set of sparse angle projection images is obtained;

[0009] Back-project the sparse angle projection map set to obtain the iterative initial value of the volume data;

[0010] Calculate the initial values ​​of the three-directional gradient fields of the volume data iteration initial value;

[0011] Based on the volume data iterative initial value and the three gradient field initial values, the forward projection and the residual back projection are calculated respectively;

[0012] Through the neural network, the data after forward projection and the data after residual back projection are regularized using three gradient field initial values;

[0013] The integral of the volume data and the gradient field is summed up to serve as the initial value for the next iteration, and the iteration is continued until completion and the result is output.

[0014] Preferably, the regularization process includes x ,D y ,D z ,Ι four directions for regularization processing; among them, D x ,D y ,D z Indicates the calculation direction gradient operator in the x, y, and z directions.

[0015] Preferably, the regularization includes performing Φ(D) expansion calculation to achieve the use of the directional gradient operator D x ,D y ,D z , expand TPO data in various directions,

[0016] The Φ(D) represents a sparse directional gradient network, D∈{D x ,D y ,D z ,I}, these four sub-networks share weights in all iterations.

[0017] Preferably, the regularization is performed by converting the image into a sparse space through an encoder, then performing feature weighting on the image in the sparse space using a shrinkage threshold function, and finally restoring the image to its original size through a decoder.

[0018] Preferably, the contraction threshold function Ι(x;μ') is expressed as: Ι(f;μ')=sign(f)max(|f|-μ',0),

[0019] Where f represents the output of the encoder, μ′ represents the contraction threshold, and the coefficient μ x , μ y , μ z Related.

[0020] Preferably, the TPO data is obtained by the following formula:

[0021] Among them, k represents the current iteration round, λ represents the iteration step size, ω k represents the k-th iteration result of ω; d refers to the gradient field in the iterative process; Ηω represents the projection data of the forward projection process in the reconstructed data; Η T ω represents back-projected projection data in the reconstructed data.

[0022] Preferably, the summing of the integral of the volume data and the gradient field refers to the summing of the integral of the volume data and the gradient field to obtain the volume data ω of the current iteration number. k :

[0023] Among them, P c (·) represents a non-negative constraint, and Denotes the directional gradient operator D in the x, y, and z directions respectively x ,D y ,D z inversion.

[0024] Preferably, before the step of adding the integrals of the volume data intermediate value and the gradient intermediate value, the method further comprises performing Nesterov acceleration processing on the three gradient fields.

[0025] Preferably, the acceleration process is calculated by the following formula:

[0026] Among them, k is the current iteration round; t k is the factor used for Nesterov acceleration; g k is a directional gradient field in x, y, or z in the current iteration round; d k =Dω k , D represents the directional gradient operator D x ,D y ,D z .

[0027] According to a second aspect of an embodiment of the present invention, a sparse atlas-based image reconstruction system is provided, comprising a processor and a memory, wherein the processor reads a computer program in the memory to implement the aforementioned sparse atlas-based image reconstruction method.

[0028] Compared with the existing technology, the image reconstruction method and system provided by the present invention can obtain the projection data required for reconstruction in a shorter time, can effectively handle the sparse angle / limited angle problem, and obtain higher-quality reconstruction results, which is particularly suitable for application in static CT systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] FIG1 is a schematic diagram of the projection structure of a static CT system;

[0030] Figure 2 is a schematic diagram of the projection operator under different projection geometries;

[0031] FIG3 is a schematic diagram of the projection angle distribution of a static CT scan;

[0032] Figure 4 is a schematic diagram of the working principle of the Directional TV Iterative Network DTV-NET;

[0033] FIG5 is a schematic diagram of the symmetrical structure of the encoder and decoder;

[0034] FIG6 is a schematic flow chart of an image reconstruction method provided by an embodiment of the present invention;

[0035] Figure 7(a) is a comparison of the effects of different image reconstruction methods under the AAPM dataset indicator evaluation;

[0036] Figure 7(b) is a comparison of the effects of different image reconstruction methods under the evaluation of the Cq500 dataset;

[0037] Figure 7(c) is a comparison of the effects of different image reconstruction methods under the evaluation of the anthropomorphic phantom index;

[0038] Figure 7(d) shows the comparison of the effects of different image reconstruction methods under the evaluation of anthropomorphic head model indicators;

[0039] FIG8 is a schematic structural diagram of an image reconstruction system provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0040] The technical content of the present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0041] As shown in FIG6 , the image reconstruction method provided by the embodiment of the present invention mainly includes the following steps.

[0042] S1: Using CT scanning, a sparse angle projection map set p is obtained.

[0043] In one embodiment of the present invention, a static CT system is used as an example for illustration. However, the image reconstruction method provided by the embodiment of the present invention is not limited to image processing obtained by scanning a static CT system, but may also be image processing of a set of sparse angle projection images obtained by other scanning methods.

[0044] In a static CT system, the radiation source is integrated into a ring-shaped structure; the detector is integrated into another ring-shaped structure. The two ring structures are arranged roughly parallel to each other, or coaxially (on the z-axis). Specifically, as shown in Figure 1, the detector ring of a static CT system is a polygonal structure composed of multiple sub-plates (close to a ring, but not a ring structure).

[0045] As we all know, projection, as a linear operation, can theoretically be constructed using a multi-layer perceptron (MLP). However, in practical applications, due to the large dimensionality of the cone-beam projection matrix and the large number of weight parameters of densely connected perceptrons, it is not feasible to apply it. To solve this problem, forward and reverse projection can be converted into a hidden layer using tensor representation, where the forward projection is the pixel drive corresponding to each detector pixel passing through the voxel and the accumulated pixel drive, while the reverse projection is the voxel drive that calculates the ratio of the corresponding pixel value to the ray length for a certain voxel. The two constitute the tensor operator of the system matrix. In the case of an ideal circular trajectory, as shown in part a of Figure 2, the ray source, the rotation center, and the center of the flat-panel detector are collinear, and the tensor operator at a certain angle can be determined by the source-to-center distance (SID) and the source-to-detector center distance (SDD). However, such existing technologies directly treat the polygonal detector structure as a ring structure for forward and reverse projection, resulting in large errors. To this end, the embodiment of the present invention adopts the forward and reverse projection model shown in part b of Figure 2, that is, a non-circular trajectory geometry with a cone angle is adopted, and the sub-detector projection tensor operator at each angle is represented by the following quantities: 1) the center point coordinate s of the detector; 2) a pair of orthogonal direction vectors passing through the center of the detector and located in the detector plane; 3) the position coordinates of each ray source.

[0046] Correspondingly, in the embodiment of the present invention, the forward and reverse projection calculations in the reconstruction algorithm are implemented as a tensorization operator (TPO), and TPO is used as a layer of a deep iterative network to complete the iterative forward and reverse projection links, thereby realizing a closed loop in the training of the reconstruction network model.

[0047] For example, using X-ray source #1 (Tube #1) as an example, when viewed in the yz plane, the X-rays emitted by it are incident on the corresponding detector's imaging area at a certain angle. In the industry, the z-axis refers to the direction in which the CT bed enters and exits the CT system; the x-axis is the horizontal axis perpendicular to the z-axis; and the y-axis is the vertical axis perpendicular to the xz plane.

[0048] Therefore, in order to adapt to the multi-source and multi-panel imaging model of the static CT system and support the oblique cone angle, it is necessary to design a specific projection tensor operator for the axial scan data of the static CT system.

[0049] There are two commonly used scanning modes for static CT systems. In the embodiments of the present invention, the first scanning mode and the second scanning mode are described respectively. Both scanning modes are based on an example in which the number of ray sources is 24 and the number of detectors is 64.

[0050] The first scanning method is a static scan with focus shift. The X-ray source and detector remain completely stationary during the scanning process. Each of the 24 X-ray sources is exposed sequentially, producing 24 projections at different angles, each forming a set of projection data. After each set of exposures is completed, the electromagnetic deflection device within the X-ray source acts on the electron beam to shift the focus of the tube for the next set of exposures. This resulting set of data is rotated slightly relative to the previous set. This process is repeated to produce 24 clusters of projection data.

[0051] The second scanning mode is a small-range rotation scan. This scanning mode is similar to the first scanning mode, except that the source and detector perform a slow mechanical rotation during the scanning process. That is, after each exposure of 24 projections, they rotate a small angle before the next exposure.

[0052] Figure 3 shows the distribution of projection angles for static CT scans. As shown in Figure 3, the circumference represents a 360° scanning angle, and a point on the circumference indicates the presence of a projection image at that angle.

[0053] Regardless of the scanning method used, each scan collects a set of projection images that are both sparse (the number of projection angles is less than 360°) and limited (the angle sampling does not evenly cover 360°). These projection image sets are denoted as p. That is, the sparse angle projection image set p has sparse angles and limited angles. The sparse angles refer to the total number of projection angles being less than 360, and the limited angles refer to the set of projection angles not evenly covering the entire 360° range.

[0054] S2: Back-project the sparse angle projection map set p to obtain the iterative initial value ω0 of the volume data.

[0055] As shown in FIG4 , by reconstructing the sparse angular projection map set p (for example, using the SIRT iterative algorithm), an iterative initial value ω0 of the volume data ω is constructed, wherein the forward projection and back projection processes in the SIRT iterative algorithm use a tensorization operator.

[0056] S3: Calculate the initial values ​​of the three-directional gradient fields of the volume data iteration initial value and

[0057] As shown in Figure 4, the three directional gradient fields and Represents the initial value of the gradient field in the x, y, and z directions (k = 0). ω0 is the initial value of the volume data, and the gradient in the x, y, and z directions is directly obtained for ω0. and

[0058] Those skilled in the art will appreciate that step S2 and step S3 may be calculated simultaneously or in reverse order, which is not a limitation and is only used as an example.

[0059] S4: Based on the iterative initial value of the volume data, the initial values ​​of the three-directional gradient fields and the sparse angle projection map set, the forward projection and the back projection are calculated to obtain the TPO data.

[0060] The TPO data calculated using the tensor quantization operator (TPO) of forward and reverse projection is shown in FIG4 , and the TPO data calculated using the tensor quantization operator (TPO) is calculated as follows:

[0061] Among them, k represents the current iteration round, λ represents the iteration step size, ω k represents the k-th iteration result of ω, Ηω represents the projection data of the forward projection process in the reconstructed data, Η T ω represents the back-projected projection data in the reconstructed data. Here, d refers to the gradient field during the iteration process, which returns to the volume data ω after the inverse operation of the directional gradient operator.

[0062] S5: Through the neural network, the TPO data is regularized in the four directions of x, y, z, and I to obtain the intermediate values ​​of the gradients in the three directions and the intermediate values ​​of the volume data.

[0063] As shown in Figure 4, the TPO data is expanded in three directions using the directional gradient operator in three directions, and the expansion is performed in the I direction to obtain the expanded data. Expansion means converting the block data itself into the original data and then adding the gradient field in three directions (i.e., expanding the data by a factor of four).

[0064] Φ(D) (which can be regarded as the Prox operator) can be understood as inputting the image into the neural network for regularization operation. After the regularization operation is performed on z, it becomes g. For example, for z k 、 Perform regularization operations (i.e., after the Φ(D) layer) to obtain the intermediate values ​​ω′ of the volume data. k and the median gradient in the three directions

[0065] Here, Φ(D) represents a network for sparsifying (or regularizing) a certain directional gradient, D∈{D x ,D y ,D z ,I}, these four sub-networks share weights in all iterations.

[0066] The extended calculation with Φ(D) means that the directional gradient operator D is calculated based on the three directions of the three-dimensional volume data. x ,D y ,D z , for the volume data Z after TPO k After expansion in all directions, further regularization of the x-direction gradient field Φ(D x ), the regularization of the y-direction gradient field Φ(D y ), the sparseness of the z-direction gradient field Φ(D z ) and its own regularization Φ(I).

[0067] It should be noted that the regularization in the embodiment of the present invention is achieved by converting the image into a sparse space through the encoder, then using a shrinkage threshold function to perform feature weighting on the image in the sparse space, and finally restoring the image to its original size through the decoder. The encoder and decoder are symmetrical. The decoder is shown in Figure 5. The encoder consists of 5 layers of 3D convolutional layers, and the size of the convolution kernel is 3×3×3. The decoder is symmetrical to the encoder and uses deconvolution to gradually restore the image. The shrinkage threshold function I(x;μ′) is the same between the encoder and decoder.

[0068] If the value of the shrinkage threshold function is less than the shrinkage threshold, it means that the voxel does not need to be thinned (regularized), and if it is greater than the shrinkage threshold, it means that the voxel needs to be thinned (regularized). These two steps can be uniformly represented by the shrinkage threshold function.

[0069] The contraction threshold function Ι(x;μ′) is expressed as: Ι(f;μ′)=sign(f)max(|f|-μ′,0) (2)

[0070] Where f represents the output of the encoder, μ′ represents the contraction threshold, and the regularization coefficient μ x , μ y , μ z Similar to the iteration step size and acceleration factor, μ′ can also be set as a learnable network parameter.

[0071] S6: Based on Nesterov acceleration processing of three gradient fields.

[0072] As shown in Figure 4, based on the three initial values ​​of the gradient field and Perform Nesterov acceleration processing (denoted by "Nes" in the figure) according to the following formula to obtain and g k+1 =g′ k +t k (g′ k -g k) (3)

[0073] Among them, k is the current iteration round; t k is the factor used for Nesterov acceleration; g k is the gradient field in one direction of x, y, or z in the current iteration (for example, when k = 0, g k for and one of the ); g′ k is the intermediate iteration variable, the intermediate value of the gradient in three directions one of the.

[0074] S7: Add the integral of the volume data intermediate value and the gradient intermediate value to serve as the initial value for the next iteration.

[0075] As shown in Figure 4, the volume data and the integral of the gradient field are summed to obtain the volume data ω of the current iteration number. k , that is, ω is calculated as follows k :

[0076] That is, calculation The harmony.

[0077] Among them, P c (·) represents a non-negative constraint, and Denotes the directional gradient operator D in the x, y, and z directions respectively x ,D y ,D z The inverse operation of .

[0078] S8: After a predetermined number of iterations, the iterated volume data is output as a result.

[0079] As can be seen, the image reconstruction method provided by the embodiments of the present invention first converts the gradient into a sparse space through an encoder, then uses a shrinkage threshold function to perform feature weighting on the image in the sparse space, and finally uses a decoder to restore the image to its original size. The encoder and decoder are symmetrical. Therefore, the image reconstruction method provided by the embodiments of the present invention achieves more accurate image pixels, better image structural information, and stronger network sparsity resistance.

[0080] For ease of understanding, the modeling process of the image reconstruction method provided by the embodiment of the present invention is further described below.

[0081] For the image reconstruction problem under static CT geometry, the image reconstruction model with direction TV as the regularization term can be used: st|D x ω|1≤δ x ,|D y ω|1≤δ y ,|D z ω|1≤δ z ,ω≥0 (6)

[0082] Where ω is the three-dimensional volume data, p represents the projection data, Η is the projection transformation matrix, D x ,D y ,D z Respectively represent the calculation of the directional gradient operator in the three directions of the three-dimensional volume data, δ x ,δ y ,δ z Represent the retention thresholds of the three directional gradients, Ηω and Η T ω represents the forward and reverse projection process in the reconstructed data, which can be represented by the tensor operator of forward and reverse projection.

[0083] For equation (6), the constraints can be eliminated using a generalized real-valued function. The Lagrangian function of the original polynomial is then calculated, and the dual problem is converted to a CP algorithm for solution. Essentially, this approach involves applying l1 regularization to the image's gradient space to find a sparse solution to the gradient. Advantageously, the present invention provides a more reasonable sparse representation of the gradient space, which facilitates better artifact removal and edge protection, while reducing the number of iterations and improving program efficiency.

[0084] use To represent the transformation space in three directions respectively, formula (6) can be expressed as:

[0085] Among them, μ x , μ y , μ z Respectively represent the sparseness of the gradient space constraints in the three directions. At the same time, based on formula (7), the one-norm transformation space of the image domain is introduced:

[0086] Where μ represents the constraint coefficient in the image space, and I represents the identity matrix. The first term in Equation (8) is the fidelity term of the conventional iterative algorithm. The second, third, and fourth terms are the regularization terms for the gradient fields in the x, y, and z directions of the volume data, while the fifth term is the regularization term for the volume data itself. This shows that in addition to the 3D volume data itself, its gradient fields in the three directions are also used as the regularization operator to improve reconstruction quality.

[0087] In the iterative reconstruction of CT with sparse angular data, sparse angular radial artifacts will inevitably be generated during the reconstruction of the volume data, which appear as divergent or windmill-shaped line artifacts. The algorithm described in the present invention can extract the characteristics of the line artifacts by performing gradient calculations on the volume data in all directions to distinguish them from the relatively smooth normal tissue structures around them. In this way, the amount of data for training the neural network can be 4 times the original amount (volume data itself + 3 gradient fields), and the artifact characteristics of the volume data are also better reflected. This is beneficial to the training of the neural network, can greatly improve the accuracy of the regularization, and thus accelerate the convergence of the iterative algorithm.

[0088] For the sake of simplicity, only the x direction is discussed here. The situations in the y and z directions are similar and will not be repeated here.

[0089] Let d = D x ω, then Formula (8) can be directly transformed into the optimization of d, and the problem can be transformed into:

[0090] For directional TV, this conversion means a faster iteration rate. satisfy According to the properties of the compact wavelet frame, equation (9) becomes:

[0091] Therefore, it is transformed into a Lasso problem model and a wavelet problem model. The FISTA algorithm is used to solve the above model and the following iterative form is obtained:

[0092] Among them, k is the current iteration round, λ is the iteration step size, t k is the factor used for Nesterov acceleration, z k ,g k are all intermediate iteration variables, P c (·) represents a non-negative constraint. Finally, the k-th iteration result of formula (11) can be expressed as ω k ,satisfy:

[0093] In order to further improve the quality of the reconstructed image, the Prox operator in Equation (11) is changed to a learnable encoder-decoder structure and extended to D x ,D y ,D z ,Ι four directions, unified as Φ(D). The iteration step λ and the momentum step t k is a fixed or learnable hyperparameter. In the figure and the following simulation analysis, the network model is called DTV-NET.

[0094] Below, the meaning of the tensor quantization operator TPO is specifically explained.

[0095] As mentioned above, ω is volume data (also called “volume block data”), and p is a set of projection image data at all angles: p = {p θ |all projection angles θ}.

[0096] For a point ω(c) in the object to be scanned under a certain projection angle θ and the corresponding detector projection point p θ (i, j), establish the equation system of its coordinates (i, j): c+t(cs)=iu s +jv s

[0097] Where t is the extension value of the ray (the above is a three-dimensional system of equations, and t, i, and j are three variables, and the values ​​of t, i, and j can be uniquely solved). Using Cramer's rule to solve the above system of equations, we get (the value of t can be solved but will not be used):

[0098] Therefore, the forward projection operator H can be expressed as the following process:

[0099] H:ω→p(θ) can be written as: θ =∫ l ω(l)dl,l=u s i+v s js

[0100] Back projection HT can be expressed as the following process

[0101] H T :{p|p(θ)forallθ}→ω can be written as:

[0102] With the forward and reverse projection operators H and HT, the entire TPO operator can be expressed as (λ is the step size): ω,p→ω+λH T (Hω-p).

[0103] The above analysis shows that a tensorization operator suitable for the stitching of multiple flat-panel detectors is proposed in the iterative forward and back projection process. The regularization term uses a directional TV deep learning model for end-to-end reconstruction of incomplete data under this static CT geometry. Experiments have demonstrated that the proposed tensorization operator effectively ensures the consistency of forward and back projections during iterative reconstruction, enabling multi-layer reconstruction with limited GPU resources. Furthermore, the DTV-NET iterative method removes sparse / limited-angle artifacts, structural artifacts, and noise, significantly improving the imaging efficiency and quality of static CT systems. The forward and back projection tensorization operators are applicable to detectors of arbitrary shapes, including polygonal detectors in static CT systems. The DTV-NET deep reconstruction algorithm is suitable for sparse and limited-angle projection data and is therefore suitable for the focal shift scanning method of static CT systems. This scanning method eliminates the need for large-scale gantry rotation and significantly improves temporal and spatial resolution.

[0104] Next, simulation data is used to illustrate the technical effects of the image reconstruction method provided by the embodiment of the present invention.

[0105] The network model DTV-NET in the embodiments of the present invention was validated using both simulation and clinical data. The simulation datasets consisted of two types: the AAPM dataset annotated by Mayo Clinics, and the Cq500 dataset. The AAPM dataset included 2,000 slices from 10 patients, with the validation set accounting for 12.7% of the total. The test volume included 14 sections, such as the lungs and abdomen. The Cq500 dataset included 1,367 slices from 23 volunteers, with the validation set accounting for 21.9% of the total. The test volume included 25 sections, such as the head and chest.

[0106] Clinical data were collected using 24-source static CT on two phantoms: an anthropomorphic human body phantom and an anthropomorphic head phantom. Each phantom collected 1080 views, with an interval of 0.33°. The structural parameters of the static CT system are as follows: R = 709 mm, r = 434.75 mm, z0 = -56.3 mm, z1 = 24.13 mm, the pixel size of each detector is 80 × 128, and the interval between each pixel is 0.265 mm. This test used the bixFDK reconstruction algorithm with a 1080-degree angle to obtain the labels (ground truth images) of the clinical data.

[0107] In the experiment, the projection angle θ={i∈24N,0≤N≤14|θ i ,θ i+1 ,θ i+2}, meaning each source is displaced 2° in axial scanning mode, resulting in three projection images. These three projection images are spaced 1° apart, forming a cluster. Therefore, the total number of projection images for the 24 sources is 72. We define views to represent the total number of projection images for the 24 sources, and intervals to represent the interval between projection angles for each source. The distribution of projection angles represented by views and intervals is shown in Figure 3. The circle represents the 360-degree scanning angle, and a point on the circle indicates the presence of a projection image at that angle.

[0108] Both the simulation and clinical datasets use the structural parameters of a static CT system. The experimental software and hardware are based on the PyTorch framework and an Nvidia A100 GPU workstation, respectively.

[0109] Figures 7(a) to 7(c) show the reconstruction results of the AAPM and Cq500 datasets and the collected phantom data using different methods. The window width of the images is [0, 1500] HU. The rMSE, SSIM, and PSNR metrics were used to evaluate the image quality of the different methods (FDK, FISTA, DTV-CP, FISTA-Net, and the DTV-NET method provided by the present invention). The results show that the image reconstruction method provided by the embodiment of the present invention (DTV-NET) achieved the highest results in all three categories. Among them, better rMSE and PSNR indicate more accurate image pixels. A higher SSIM indicates better image structural information and stronger network anti-sparseness capabilities.

[0110] Based on the above-mentioned sparse atlas-based image reconstruction method, the present invention also provides a sparse atlas-based image reconstruction system. This image reconstruction system can be a static CT system, or it can be security inspection equipment, non-destructive testing equipment, etc. As shown in Figure 8, the image reconstruction system includes one or more processors 21 and at least one memory 22. The memory 22 is coupled to the processor 21 and is used to store one or more programs. When the one or more programs are executed by the one or more processors 21, the one or more processors 21 implement the sparse atlas-based image reconstruction method as described in the above-mentioned embodiment.

[0111] The processor 21 is used to control the overall operation of the image reconstruction system to complete all or part of the steps of the above-mentioned sparse atlas-based image reconstruction method. The processor 21 can be a central processing unit (CPU), a graphics processing unit (GPU), a field programmable gate array (FPGA), an application-specific integrated circuit (ASIC), a digital signal processing (DSP) chip, etc. The memory 22 is used to store various types of data to support the operation of the image reconstruction system. Such data may include, for example, instructions for any application or method operating on the image reconstruction system, as well as application-related data. The memory 22 can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic memory, flash memory, etc.

[0112] In an exemplary embodiment, the image reconstruction system can be specifically implemented by a computer chip or entity, or by a product with a certain function, for executing the above-mentioned sparse atlas-based image reconstruction method and achieving the same technical effect as the above-mentioned method.

[0113] In another exemplary embodiment, the present invention further provides a computer-readable storage medium comprising program instructions, which, when executed by a processor, implements the steps of the sparse atlas-based image reconstruction method in any one of the above embodiments.

[0114] The order of the steps in the above embodiments can be adjusted according to actual needs, and other steps can be inserted or added, such as pre-processing the volume data. Moreover, the calculation formula can be replaced by other formulas as long as the technical purpose of each step is met.

[0115] The above describes in detail the sparse atlas-based image reconstruction method and system provided by the present invention. For those skilled in the art, any obvious modification to the present invention without departing from its essence would constitute an infringement of the present invention's patent rights and would incur corresponding legal liability.

Claims

1. A sparse atlas-based image reconstruction method, characterized in that The following steps are involved: Using CT scanning, a set of sparse angle projection images is obtained; Back-project the sparse angle projection map set p to obtain the iterative initial value of the volume data; Calculate the initial values ​​of the gradient fields in three directions for the initial values ​​of the volume data iteration; Based on the volume data iteration initial value, the three-directional gradient field initial value and the sparse angle projection map set, the forward projection and the reverse projection are calculated to obtain the TPO data; Through the neural network, the TPO data is regularized in four directions of x, y, z, and Ι to obtain the intermediate values ​​of the gradients in the three directions and the intermediate values ​​of the volume data, where x, y, and z represent the three directions in the three-dimensional space and Ι represents the unit matrix; The integral of the volume data intermediate value and the gradient intermediate value is summed up to serve as the initial value for the next iteration, and the iteration is repeated until completion and the result is output.

2. The image reconstruction method based on sparse atlas according to claim 1, characterized in that: The regularization process includes x ,D y ,D z , Ι regularization is performed in four directions, Where D x ,D y ,D z Represents the calculation direction gradient operator in the x, y, and z directions.

3. The image reconstruction method based on sparse atlas according to claim 2, characterized in that: The regularization includes performing a Φ(D) expansion calculation to achieve the use of the directional gradient operator D x ,D y ,D z , expand TPO data in all directions, The Φ(D) represents a sparse directional gradient network, D∈{D x ,D y ,D z ,I}, these four sub-networks share weights in all iterations.

4. The image reconstruction method based on sparse atlas according to any one of claims 1 to 3, characterized in that: The regularization is performed by converting the image into a sparse space through an encoder, then using a shrinkage threshold function to perform feature weighting on the image in the sparse space, and finally restoring the image to its original size through a decoder.

5. The image reconstruction method based on sparse atlas according to claim 4, characterized in that: The contraction threshold function Ι(x;μ') is expressed as: Ι(f;μ')=sign(f)max(|f|-μ',0), Where f represents the output of the encoder, μ′ represents the contraction threshold, and the coefficient μ x , μ y , μ z Related.

6. The image reconstruction method based on sparse atlas according to any one of claims 1 to 3, characterized in that: The TPO data is obtained by the following formula: With k =ω k +λH T (Hω k -p) Among them, k represents the current iteration round, λ represents the iteration step size, ω k represents the k-th iteration result of ω; d refers to the gradient field in the iterative process; Ηω represents the projection data of the forward projection process in the reconstructed data; Η T ω represents back-projected projection data in the reconstructed data.

7. The image reconstruction method based on sparse atlas according to claim 6, characterized in that: The summing up of the integral of the volume data and the gradient field refers to summing up the integral of the volume data and the gradient field to obtain the volume data ω of the current iteration number. k : Among them, P c (·) represents a non-negative constraint, and Denote the directional gradient operator D in the x, y, and z directions respectively x ,D y ,D z inversion.

8. The image reconstruction method based on sparse atlas according to claim 7, characterized in that Before the step of adding the integrals of the volume data intermediate value and the gradient intermediate value, the method further includes performing Nesterov acceleration processing on the three gradient fields.

9. The image reconstruction method based on sparse atlas as claimed in claim 8, characterized in that The acceleration process is calculated by the following formula Where k is the current iteration round; t k is the factor used for Nesterov acceleration; g k is a directional gradient field in x, y, or z of the current iteration round; d k =Dω k , D represents the directional gradient operator D x ,D y ,D z .

10. An image reconstruction system based on sparse atlas, characterized in that It comprises a processor and a memory, wherein the processor reads a computer program in the memory to implement the image reconstruction method based on a sparse atlas as claimed in any one of claims 1 to 9.

Images