Sparse view angle 3D-DSA reconstruction method based on three-dimensional Poisson generative model

Abstract

The invention discloses a sparse view angle 3D-DSA reconstruction method based on a three-dimensional Poisson generative model. The method comprises the following steps: obtaining pairing data of a sparse view angle 2D-DSA projection drawing and a 3D-DSA reconstruction image; extracting features of a projection image by using a projection domain encoder, then converting features of two-dimensional projection into a three-dimensional image domain according to a geometrical relationship of the cone beam CT, and then obtaining a prior image by using an image decoder; the method comprises the following steps: constructing a three-dimensional Poisson generation model, adding noise to a three-dimensional image patch in a training stage to obtain a disturbance image, outputting a network reconstruction image by taking a prior image as a condition and the three-dimensional image patch before noise addition as a target image, and calculating loss of the output image and the target image to update network parameters; the mean square error loss and the mean square error loss of the maximum intensity projection images of the three orthogonal planes are used during loss calculation; in the sampling stage, random noise is used as input, a prior image is used as a condition, noise of a noise image is continuously denoised within a limited step length, and finally a reconstructed 3D-DSA image is obtained.

Description

Technical Field

[0001] The present invention belongs to the technical field of medical image processing, and relates to a sparse view medical image reconstruction technology, and in particular to a sparse view 3D-DSA reconstruction method based on a three-dimensional Poisson generation model. Background Art

[0002] Cardiovascular and cerebrovascular diseases are common and frequently occurring diseases and are one of the leading causes of death. Seeking more effective diagnostic and treatment methods for cardiovascular and cerebrovascular diseases is of great practical significance in reducing their morbidity and mortality, and alleviating the pressure on the healthcare system. Digital subtraction angiography (DSA) is an interventional surgical technique that uses X-ray imaging technology to image blood vessels. It can provide high-quality and detailed vascular images and is widely used in interventional clinical practice. It is considered the gold standard for the diagnosis of acute cerebrovascular diseases.

[0003] Conventional anteroposterior and lateral DSA are often affected by the overlapping angles of multiple vessels and improper imaging angle selection, blurring the relationship between the diseased and surrounding vessels and even misdiagnosing vascular stenosis due to occlusion. To address these issues, three-dimensional rotational digital subtraction angiography (3D-DSA) has emerged. 3D-DSA can visualize vascular images from multiple angles with high accuracy, clearly demonstrating vascular stenosis and possessing significant clinical value. However, in current commercial DSA systems, 3D-DSA reconstruction still relies on the traditional FDK (filtered back projection) algorithm, typically requiring hundreds of scan views to produce artifact-free, high-quality images. This large number of X-ray projections results in significant radiation exposure for both patients and radiologists. Reducing the number of scans results in reconstructions containing many artifacts and broken vascular structures. Therefore, reconstructing high-quality 3D-DSA images using a small number of projections is crucial.

[0004] The current 3D-DSA reconstruction algorithm is based on the cone-beam CT reconstruction algorithm. Therefore, generally speaking, the research results of the cone-beam CT reconstruction algorithm can be directly applied to the 3D-DSA reconstruction algorithm. The current cone-beam CT reconstruction algorithms can be mainly divided into the following categories: iterative reconstruction algorithms and deep learning-based reconstruction algorithms. The iterative reconstruction method has the following problems: first, the performance of the iterative algorithm is heavily dependent on the design of the regularization term and the selection of hyperparameters, which usually require a lot of experience and experiments to adjust; second, the iterative reconstruction method has a high computational cost. Although it does not rely on GPU resources, it has certain requirements for CPU resources; third, although its reconstruction results are less noisy than those of FDK, vascular details are also lost.

[0005] Reconstruction algorithms based on deep learning mainly include image domain methods, dual-domain methods, diffusion model-based methods, and differentiable rendering-based methods. (1) Image domain methods usually use the FDK algorithm to obtain an initial result and then use a neural network to restore it. In this process, artifacts and noise are introduced into the reconstructed image, affecting the ability of subsequent algorithms to distinguish structural features. On the other hand, since the cross-sectional data of 3D-DSA in the image domain often presents point-like and polygonal block features, the blood vessels in the cross section do not have a significant regular structure under the influence of random noise and streak artifacts. (2) Among the strategies of dual-domain methods, a widely used approach is to use two cascaded networks to simultaneously process data in the projection domain and image domain. However, these methods often require all two-dimensional projection sequences to be input in a single inference, which will take up a lot of computing resources in network training and practical applications. Due to this limitation, such methods often use MLP to fuse feature data from different projection domains in practical applications. On the one hand, this design only utilizes the projection domain features and abandons the cross-utilization of the image domain structural features, so the reconstruction quality is still limited. On the other hand, such methods can only process projection views with low resolution, and are therefore limited in high-resolution DSA images. (3) Generative models based on physics inspiration have achieved remarkable results in CT imaging tasks, especially the diffusion model and Poisson flow generation model, which have shown good performance in unconditional and conditional image generation. Both types of generative models are based on physics inspiration: the core idea of ​​the diffusion model is derived from non-equilibrium thermodynamics, while the Poisson flow generation model is derived from electrostatics. However, most methods based on the diffusion model are trained from two-dimensional slice images, which is not conducive to the structure of three-dimensional blood vessels; the sampling step size is large and time-consuming; and the methods based on the Poisson flow diffusion model are currently mainly used in two-dimensional image denoising. (4) Differentiable methods mainly include NERF and Gaussian sputtering. However, these methods require individual training for each patient. Using NERF to reconstruct a single CBCT image takes tens of minutes or even two hours, and they do not completely overcome the challenges brought by extremely sparse views. Although the Gaussian sputtering method improves the reconstruction speed, it still takes a long time to obtain a high-quality reconstructed image due to the sparsity of vascular information.

[0006] In summary, how to combine the dual-domain information of high-resolution projection domain 2D-DSA images and image domain 3D-DSA images to obtain high-quality reconstruction results in the case of sparse view angles is still a problem that needs to be studied. Summary of the Invention

[0007] Objective of the Invention: This invention aims to provide a sparse-view 3D-DSA reconstruction method based on a 3D Poisson generative model. First, a projection information extraction model is constructed to obtain prior conditions from projection images. The model consists of three components: a projection domain image encoder, a feature transformation, and an image domain decoder. The projection domain encoder extracts features from the input N projection domains, then transforms the features to the image domain based on the geometric relationships of cone-beam CT. Finally, the image domain decoder converts the features into an image. The projection information extraction model is trained on z-axis slices of the image, using the L1 loss to update parameters. After obtaining the prior image, a 3D conditional Poisson generative model is constructed to maintain the continuity of the 3D vascular structure. The model is trained and sampled in the 3D image space. During the model training phase, the input is a noise-perturbed image patch and a prior image patch. The output is the reconstructed image patch. The parameters are updated by computing the mean square error loss and the maximum intensity projection loss on three orthogonal planes with the reference image. During the model sampling phase, the input is random noise and the prior image patch. At each step, the current denoised image is output and corrected using the Euler method to obtain the final 3D-DSA reconstruction result.

[0008] To achieve the above object, the present invention provides a sparse view 3D-DSA reconstruction method based on a three-dimensional Poisson generation model, comprising the following steps:

[0009] S1: Obtain N sparse-view 2D projection images and 3D-DSA image I reconstructed using the FDK algorithm based on full-view projection as the gold standard;

[0010] S2: For each z-axis slice of the 3D-DSA image, the projection domain image range [0, h max ;0,w max ];

[0011] S3: Construct a projection information extraction model. For each z-axis slice of the 3D-DSA image, first use the projection domain image encoder E to extract the features f of N two-dimensional projection images according to the image range obtained in step S2. proj , and then the features of the projection image are converted into image domain features f according to the geometric relationship of cone-beam CT slice ; Then construct the image domain decoder D, and transform the image domain feature f slice Decode into image slice';

[0012] S4: Construct an L1 loss function based on the slice in step S2 and the slice' in step S3, and use the optimizer to update the parameters of the projection information extraction model;

[0013] S5: Obtain a priori image C based on the projection information extraction model trained in step S4 and the two-dimensional projection image of step S1; perform image segmentation on the priori image C and the gold standard image I to obtain a priori image block C patch and gold standard image patch I patch ;

[0014] S6: Construct a three-dimensional conditional Poisson generative model. During the model training phase, the gold standard image block I obtained in step S5 is patch Noise is added, and the model input is the noisy perturbation image block I patch’ and the prior image C as the guiding condition patch’ , the model output is the reconstructed image block I net ;

[0015] S7: According to steps S5 and S6, patch , I net Construct the mean square error loss function L mse and the maximum intensity projection loss function L on three orthogonal planes proj ,constraining the consistency between the model output image and the gold standard image at different observation angles,improving the accuracy of the results;

[0016] S8: The mean square error loss function L in step S7 is mse With the maximum intensity projection loss L proj Perform weighted calculations and use an optimizer to update parameters of the three-dimensional conditional Poisson generation model;

[0017] S9: In the model sampling stage, the three-dimensional conditional Poisson generative model trained in step S8 and the prior image C in step S5 are used. patch , with random noise and C patch As the input of the model, the noise is continuously denoised within a limited step size to obtain the final reconstructed image I net ; will I net The final 3D-DSA reconstructed image is obtained by stitching.

[0018] Furthermore, the process of calculating the projection domain image range in step S2 includes constructing a projection matrix, image position point mapping, and maximum value calculation. For the point image_index(x, y, z) on the current slice image, the rotation angle of the projection image is angle, and the voxel spacing is v s , the projection pixel spacing is p s , the distance from the focus to the rotation center is SOD, the distance from the focus to the detector is SID, and the projection matrix M is constructed as

[0019]

[0020] Where cos and sin represent the cosine function and sine function, then the scaling coordinate point of the current angle is image_index×v s ×M T , the magnification coeff is SID / (SOD-z), and the final corresponding projection domain coordinates are point×coeff / p s A projection operation is performed on each point of the current slice image, and the maximum value of the projection domain coordinates is calculated. This is used to calculate the range of the projection image that actually participates in each training. This step can solve the problem of high computing resources occupied when extracting high-resolution projection image features, thereby fully extracting the features of the projection image without sacrificing image resolution.

[0021] Furthermore, the projection domain image encoder E in step S3 uses U-Net as the backbone network. U-Net is a classic medical image segmentation network. The use of U-Net can ensure the model's ability to extract medical image features. The projection domain image encoder E includes 4 downsampling modules consisting of convolutional layers, Relu activation function layers, and maximum pooling layers, and 4 upsampling modules consisting of convolutional layers and transposed convolutional layers. The input is the processed projection image of N angles, and the output is the corresponding projection image feature f proj ; The feature conversion process in step 3 includes constructing a projection matrix, image position point mapping and image interpolation process, wherein the operations of constructing the projection matrix and the image position point are consistent with those in step S2, and the image interpolation operation is because the result of projecting the points in the image domain to the projection domain is a decimal, and bilinear interpolation is required to obtain the final feature value; the image domain decoder D in step 3 also uses U-Net as the backbone network, and also includes 4 downsampling modules consisting of convolutional layers, Relu activation function layers and maximum pooling layers, and 4 upsampling modules consisting of convolutional layers and transposed convolutional layers; the input is the converted image feature f slice , the output is the decoded image slice'.

[0022] Furthermore, the L1 loss function expression constructed in step S4 is:

[0023]

[0024] Furthermore, the process of obtaining the prior image C in step S5 is to input N two-dimensional projection images into the model trained in step S4, and then stack each output slice image on the z-axis; the image cutting process in step 5 is to traverse the coordinates of the corresponding prior image and gold standard image, and cut them with a step size of 64×64×64 and a size of 128×128×128 to obtain the prior image block and the gold standard image block.

[0025] Furthermore, the three-dimensional conditional Poisson generative model in step S6 uses 3D-UNet as the backbone network, and the training phase of the model includes:

[0026] (1) Sample σ from the data distribution p(σ), calculate r = σ√D (D is the augmented dimension of the Poisson generative model), and then sample the Beta distribution S with shape parameters N / 2 (N is the size of the image block) and D / 2 beta ~Beta(N / 2,D / 2), then we get the inverse Beta distribution S inv-beta =S beta / (1-S beta ), combining the radius r and the inverse Beta distribution sample to obtain the noise perturbation radius R = r × S inv-beta Then, we sample from the standard Gaussian distribution and normalize it to obtain the direction vector of the noise perturbation v = u / ||u||, where u ~ N(0, I), and N represents the Gaussian distribution. Finally, we multiply the noise perturbation radius and the direction vector to obtain the final perturbation noise Rv.

[0027] (2) In each iteration, image block I is taken from the gold standard image. patch , add random noise to disturb and get the noisy image I patch’ =I patch +Rv;

[0028] (3) The image block I after adding noise patch’ and the corresponding prior image block C patch As the input of the 3D conditional Poisson generative model, the output is the reconstructed image block I net .

[0029] This step constructs a three-dimensional Poisson generative model and learns the inter-domain mapping from defective structures to complete structures in the three-dimensional image domain, thereby optimizing the inter-layer discontinuity problem caused by reconstruction in two-dimensional slices and improving the quality of the final reconstructed image.

[0030] Furthermore, the maximum intensity projection loss L in step S7 is proj The calculation method is to perform maximum intensity projection on the reconstructed image block and the gold standard image block in the x, y, and z planes, and obtain the maximum projection image I of the reconstructed image block in the x-axis direction respectively. net_x , reconstruct the maximum projection image I of the image block in the y-axis direction net_y , reconstruct the maximum projection image I of the image block in the z-axis direction net_z and the maximum projection image I of the gold standard image block in the x-axis direction patch_x , the maximum projection image I of the gold standard image block in the y-axis direction patch_y , the maximum projection image I of the gold standard image block in the z-axis direction patch_z, then calculate the mean square error (MSE) loss between each pair of them, and the expression is as follows:

[0031]

[0032] This step can constrain the consistency between the model output image and the gold standard image at different observation angles, thereby improving the quality of the reconstructed image.

[0033] Furthermore, during the model training in step S8, the mean square error loss and the maximum intensity projection loss are calculated using the denoised image patch and the gold standard image patch output by the model to construct a loss function. The weighted loss expression is as follows, where λ is the weighting coefficient and takes 0.1:

[0034] L total =L mse +λL proj

[0035] Furthermore, the sampling stage of the three-dimensional conditional Poisson generation model in step S9 includes:

[0036] (1) Set the maximum scale σ of the noise max , calculate r = σ√D (D is the augmented dimension of the Poisson generation model), and then sample the Beta distribution S with shape parameters N / 2 (N is the size of the image patch) and D / 2 beta ~Beta(N / 2,D / 2), then obtain the inverse Beta distribution S inv-beta =S beta / (1 - S beta ), combine the radius r and the inverse Beta distribution sample to obtain the noise perturbation radius R = r×S inv-beta ; then sample from the standard Gaussian distribution and normalize to obtain the direction vector v of the noise perturbation, where u~N(0,I), and N represents the Gaussian distribution. Finally, multiply the noise perturbation radius and the direction vector to obtain the initial noise x 0 =Rv;

[0037] (2) Set the sampling step size T = 8;

[0038] (3) For the current step size n (0 < n < T), regard the sampling process of the model as solving an ordinary differential equation using the Euler method. First, according to the current image x n , the prior image C, the noise scale t corresponding to the current step n , and the trained model f to obtain the instantaneous change rate d n :

[0039]

[0040] Then, according to the noise scale t + 1 corresponding to the next stepn and the current image x n and the instantaneous rate of change d n Get the next step length image x n+1 :

[0041] x n+1 =x n +(t n+1 -t n )d n

[0042] Then calculate the instantaneous rate of change d at t+1 ’n And correct the results:

[0043]

[0044] (4) The image blocks of each sampling result are spliced ​​according to the spatial position to obtain the 3D-DSA reconstructed image.

[0045] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0046] 1. This paper comprehensively utilizes projection domain information and image domain information and proposes a sparse view 3D-DSA reconstruction method based on a three-dimensional Poisson generative model. When a sparse view projection image is input, a high-quality 3D-DSA reconstructed image is obtained.

[0047] 2. This invention first trains a projection information extraction model to fully utilize the information of images projected at different angles and improve subsequent reconstruction results. To address the problem of high computational resources required to extract features from high-resolution projection images, the training process uses image slices as units to calculate the range of feature images to be extracted each time. Furthermore, the gradients of the projected feature images are grouped and returned in each iteration, effectively extracting features from the projection images without sacrificing image resolution and reducing computational resources.

[0048] 3. Based on the projection information extraction model, the present invention constructs a 3D Poisson generative model. This model can learn the inter-domain mapping from defective structures to complete structures in the 3D image domain, thereby optimizing the inter-layer discontinuity problem caused by reconstruction based on 2D slices and improving the quality of the reconstructed image.

[0049] 4. The three-dimensional Poisson generation model proposed in the present invention adopts the mean square error loss and the maximum intensity projection loss of three orthogonal planes, thereby constraining the consistency of the reconstructed image and the real image in different planes, thereby improving the reconstruction quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1A schematic diagram of the framework topology structure of sparse perspective 3D-DSA reconstruction based on a three-dimensional Poisson generation model provided by the present invention.

[0051] Figure 2 This is a schematic diagram of the topological structure of the projection information extraction model provided by the present invention.

[0052] Figure 3 A schematic flow chart of a sparse perspective 3D-DSA reconstruction method based on a three-dimensional Poisson generation model provided by the present invention.

[0053] Figure 4 a—j are comparison results of the reconstruction results of the test image on the simulated projection dataset.

[0054] Figure 5 a—j are comparison results of the reconstruction results of the test image on the real projection dataset. DETAILED DESCRIPTION

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

[0056] Example: In current commercial DSA systems, 3D-DSA reconstruction still relies on the traditional FDK (filtered back projection) algorithm, which typically requires hundreds of scan views to produce artifact-free, high-quality images. This large number of X-ray projections results in significant radiation exposure for both patients and radiologists. Current reconstruction methods, however, cannot effectively extract features from high-resolution vascular projection images and typically perform reconstruction in units of two-dimensional slices, which is detrimental to maintaining the continuity of vascular structure and restoring detailed details. Therefore, this method first constructs a projection information extraction model to obtain prior conditions from projection images. The model consists of three components: a projection domain encoder, a feature transform, and an image domain decoder. The projection domain encoder extracts features from the input N projection domains, then transforms the features into the image domain based on the geometric relationships of cone-beam CT. Finally, the image domain decoder converts the features into an image. The projection information extraction model is trained in units of z-axis slices, using an L1 loss to update parameters. After obtaining the prior image, a three-dimensional conditional Poisson generative model is constructed to maintain the continuity of the three-dimensional vascular structure. The model is trained and sampled in the three-dimensional image space. During the model's training phase, the input is a noise-perturbed image patch and a prior image patch. The output is a reconstructed image patch. The reconstructed image patch is compared with the gold standard image and the mean square error loss and maximum intensity projection loss on three orthogonal planes are calculated to update the parameters. During the model's sampling phase, the input is random noise and the prior image patch. At each step, the current denoised image is output and corrected using the Euler method to obtain the final 3D-DSA reconstruction result.

[0057] like Figure 3 FIG. 1 is a flow chart of the present invention, which shows a sparse view 3D-DSA reconstruction method based on a three-dimensional Poisson generation model, comprising the following steps:

[0058] S1: Obtain N sparse-view 2D projection images and 3D-DSA image I reconstructed using the FDK algorithm based on full-view projection as the gold standard;

[0059] S2: For each z-axis slice of the 3D-DSA image, the projection domain image range [0, h max ;0,w max ];

[0060] The process of calculating the projection domain image range includes constructing the projection matrix, image position point mapping and maximum value calculation. For the point image_index(x,y,z) on the current slice image, the rotation angle of the projection image is angle and the voxel spacing is v s , the projection pixel spacing is p s , the distance from the focus to the rotation center is SOD, the distance from the focus to the detector is SID, and the projection matrix M is constructed as

[0061]

[0062] Where cos and sin represent the cosine function and sine function, then the scaling coordinate point of the current angle is image_index×v s ×M T , the amplification factor coeff is SID / (SOD-z),

[0063] The final corresponding projection domain coordinates (h, w) are point×coeff / p s This operation is performed for each point of the current slice image to calculate the maximum value of the projection domain coordinates, which is used to calculate the range of the projection image that actually participates in training each time, thereby reducing the computing resource consumption during training.

[0064] S3: Construct a projection information extraction model. For each slice of the 3D-DSA image, first use the projection domain image encoder E to extract the features f of N two-dimensional projection images according to the image range obtained in step S2. proj , and then the features of the projection image are converted into image domain features f according to the geometric relationship of cone-beam CT slice ; Then construct the image domain decoder D, and transform the image domain feature f slice Decode into image slice';

[0065] The projection domain image encoder E uses U-Net as the backbone network. U-Net is a classic medical image segmentation network. Using U-Net can ensure the model's ability to extract medical image features. The encoder E contains four downsampling modules consisting of convolutional layers, ReLU activation function layers, and maximum pooling layers, and four upsampling modules consisting of convolutional layers and transposed convolutional layers. The input is the processed projection image at N angles, and the output is the corresponding projection image feature f proj The feature conversion process includes constructing a projection matrix, image position point mapping and image interpolation process, wherein the operations of constructing the projection matrix and the image position point are consistent with the step S2, and the image interpolation operation is because the result of projecting the points in the image domain to the projection domain is a decimal, and bilinear interpolation is required to obtain the final feature value; the image domain decoder D also uses U-Net as the backbone network, and also includes 4 downsampling modules consisting of convolutional layers, Relu activation function layers and maximum pooling layers, and 4 upsampling modules consisting of convolutional layers and transposed convolutional layers; the input is the converted image feature f slice , the output is the decoded image slice';

[0066] S4: Construct an L1 loss function based on the slice in step S2 and the slice' in step S3, and use the optimizer to update the parameters of the projection information extraction model;

[0067]

[0068] S5: Based on the projection information extraction model trained in step S4 and the two-dimensional projection image of step S1, a priori image C is obtained; the priori image C and the gold standard image I are cut into blocks to obtain C patch and I patch ;

[0069] The process of obtaining the prior image C is to input N two-dimensional projection images into the model trained in step S4, and then stack each output slice image on the z-axis; the image segmentation process in step 5 is to traverse the coordinates of the corresponding prior image and gold standard image, and crop them with a step size of 64×64×64 and a size of 128×128×128 to obtain the prior image block and the gold standard image block;

[0070] S6: Construct a three-dimensional conditional Poisson generative model. During the model training phase, the gold standard image block I obtained in step S5 is patch Noise is added, and the model input is the noisy perturbation image block I patch’ and the prior image C as the guiding condition patch’ , the model output is the reconstructed image block I net ;

[0071] The 3D conditional Poisson generative model uses 3D-UNet as the backbone network. The training phase of the model includes:

[0072] (1) Sample σ from the data distribution p(σ), calculate r = σ√D (D is the augmented dimension of the Poisson generative model), and then sample the Beta distribution S with shape parameters N / 2 (N is the size of the image block) and D / 2 beta ~Beta(N / 2,D / 2), then we get the inverse Beta distribution S inv-beta =S beta / (1-S beta ), combining the radius r and the inverse Beta distribution sample to obtain the noise perturbation radius R = r × S inv-beta Then, we sample from the standard Gaussian distribution and normalize it to obtain the direction vector of the noise perturbation v = u / ||u||, where u ~ N(0, I), and N represents the Gaussian distribution. Finally, we multiply the noise perturbation radius and the direction vector to obtain the final perturbation noise Rv.

[0073] (2) In each iteration, image block I is taken from the gold standard image. patch, add random noise to disturb and get the noisy image I patch’ =I patch +Rv;

[0074] (3) The image block I after adding noise patch’ and the corresponding prior image block C patch As the input of the 3D conditional Poisson generative model, the output is the reconstructed image block I net ;

[0075] S7: According to steps S5 and S6, patch , I net Construct the mean square error loss function L mse and the maximum intensity projection loss function L on three orthogonal planes proj ,constraining the consistency between the model output image and the gold standard image at different observation angles,improving the accuracy of the results;

[0076] Maximum intensity projection loss L proj The calculation method is to perform maximum intensity projection on the reconstructed image block and the gold standard image block in the x, y, and z planes to obtain the maximum projection image I of the reconstructed image block in the x-axis direction. net_x , reconstruct the maximum projection image I of the image block in the y-axis direction net_y , reconstruct the maximum projection image I of the image block in the z-axis direction net_z and the maximum projection image I of the gold standard image block in the x-axis direction patch_x , the maximum projection image I of the gold standard image block in the y-axis direction patch_y , the maximum projection image I of the gold standard image block in the z-axis direction patch_z , and then calculate the mean square error (MSE) loss between them, thereby constraining the consistency of the model output image and the gold standard image at different observation angles and improving the accuracy of the results. The calculation expression is as follows:

[0077]

[0078] S8: The mean square error loss function L in step S7 is mse With the maximum intensity projection loss L proj Perform weighted calculations and use the optimizer to update the parameters of the three-dimensional conditional Poisson generation model. The weighted loss expression is as follows, where λ is the weighting coefficient, which is 0.1;

[0079] L total =L mse +λL proj

[0080] S9: In the model sampling stage, the three-dimensional conditional Poisson generative model trained in step S8 and the prior image block C in step S5 are used.patch , with random noise and C patch as the input to the model, continuously denoise the noise within a finite number of steps to obtain the final reconstructed image patch I net ; Stitch I net to obtain the final 3D-DSA reconstructed image.

[0081] The sampling stage of the three-dimensional conditional Poisson generation model includes:

[0082] (1) Set the maximum scale σ of the noise max , calculate r = σ√D (D is the augmented dimension of the Poisson generation model), then sample the Beta distribution S with shape parameters N / 2 (N is the size of the image patch) and D / 2 beta ~Beta(N / 2,D / 2), then obtain the inverse Beta distribution S inv-beta = S beta / (1 - S beta ), combine the radius r and the inverse Beta distribution sample to obtain the noise perturbation radius R = r×S inv-beta ; Then sample and normalize from the standard Gaussian distribution to obtain the direction vector v of the noise perturbation, v = u / ||u||, where u~N(0,I), N represents the Gaussian distribution. Finally, multiply the noise perturbation radius and the direction vector to obtain the initial noise x 0 = Rv;

[0083] (2) Set the sampling step size T = 8;

[0084] (3) For the current step size n (0 < n < T), regard the sampling process of the model as solving an ordinary differential equation using the Euler method. First, according to the current image x n , the prior image C, the noise scale t corresponding to the current step size n and the trained model f to obtain the instantaneous change rate d n :

[0085]

[0086] Then, according to the noise scale t n+1 corresponding to the next step size and the current image x n and the instantaneous change rate d n to obtain the image x n+1 of the next step size:

[0087] x n+1 = x n +(t n+1 -t n )d n

[0088] Then calculate the instantaneous change rate d at t + 1’n And correct the results:

[0089]

[0090] (4) The image blocks of each sampling result are spliced ​​according to the spatial position to obtain the 3D-DSA reconstructed image.

[0091] In order to prove the effectiveness of the present invention, the present invention also provides the following experiments:

[0092] Specifically, we conducted experiments using two datasets: simulated projection and real projection. The simulated projection dataset contained 44 data points, with a training / validation / test split of 34 / 5 / 5; the simulated projection dataset contained 18 data points, with a training / validation / test split of 12 / 3 / 3. The model was trained and tested using 20 and 30 input viewpoints, respectively.

[0093] In comparative experiments, the reconstruction results of our method were compared with those of FDK, OSSART, FBP-ConvNet, DIF-Net, NAF, 3D Gaussian (3DGS), GAN, and DDPM on simulated and real projection datasets, using PSNR, SSIM, and LPIPS as evaluation metrics. As shown in Table 1, our method significantly outperforms other domain adaptation methods in terms of reconstruction performance evaluation metrics.

[0094]

[0095]

[0096] Table 1 Comparison of the reconstruction results of the sparse view 3D-DSA reconstruction method based on the 3D Poisson generative model and other reconstruction methods on simulated projection data and real projection data

[0097] In order to more intuitively demonstrate the effectiveness of this method, the visual reconstruction results on the simulated projection dataset are shown as follows: Figure 4 As shown in (aj). Where:

[0098] (a) is the true reconstruction result of the original image;

[0099] (b) is the reconstruction result of the present invention;

[0100] (c) is the reconstruction result of the FDK algorithm;

[0101] (d) is the reconstruction result of the OSSART algorithm;

[0102] (e) is the reconstruction result of the FBP-ConvNet algorithm;

[0103] (f) is the reconstruction result of the DIF-Net algorithm;

[0104] (g) is the reconstruction result of the NAF algorithm;

[0105] (h) is the reconstruction result of 3D Gaussian algorithm;

[0106] (i) is the reconstruction result of the GAN algorithm;

[0107] (j) is the reconstruction result of the DDPM algorithm;

[0108] The visual reconstruction results on the real projection dataset are shown in the figure Figure 5 As shown in (aj). Where:

[0109] (a) is the true reconstruction result of the original image;

[0110] (b) is the reconstruction result of the present invention;

[0111] (c) is the reconstruction result of the FDK algorithm;

[0112] (d) is the reconstruction result of the OSSART algorithm;

[0113] (e) is the reconstruction result of the FBP-ConvNet algorithm;

[0114] (f) is the reconstruction result of the DIF-Net algorithm;

[0115] (g) is the reconstruction result of the NAF algorithm;

[0116] (h) is the reconstruction result of 3D Gaussian algorithm;

[0117] (i) is the reconstruction result of the GAN algorithm;

[0118] (j) is the reconstruction result of the DDPM algorithm;

[0119] It can be seen that compared with other reconstruction algorithms, the method proposed in the present invention can reconstruct 3D-DSA images more accurately, restore better details, and have high-quality visualization results.

[0120] It should be noted that the above embodiments are not intended to limit the scope of protection of the present invention, and equivalent changes or substitutions made on the basis of the above technical solutions fall within the scope of protection of the claims of the present invention.

Claims

1. A sparse view 3D-DSA reconstruction method based on a three-dimensional Poisson generation model, characterized in that: The method comprises the following steps: S1: Obtain N sparse-view 2D projection images and 3D-DSA image I reconstructed using the FDK algorithm based on full-view projection as the gold standard; S2: For each z-axis slice of the 3D-DSA image, the projection domain image range [0, h max ;0,w max ]; S3: Construct a projection information extraction model. For each z-axis slice of the 3D-DSA image, first use the projection domain image encoder E to extract the features f of N two-dimensional projection images according to the image range obtained in step S2. proj , and then the features of the projection image are converted into image domain features f according to the geometric relationship of cone-beam CT slice ; Then construct the image domain decoder D, and transform the image domain feature f slice Decode into image slice'; S4: Construct an L1 loss function based on the slice in step S2 and the slice' in step S3, and use the optimizer to update the parameters of the projection information extraction model; S5: Obtain a priori image C based on the projection information extraction model trained in step S4 and the two-dimensional projection image of step S1; perform image segmentation on the priori image C and the gold standard image I to obtain a priori image block C patch and gold standard image patch I patch ; S6: Construct a three-dimensional conditional Poisson generative model. During the model training phase, the gold standard image block I obtained in step S5 is patch Noise is added, and the model input is the noisy perturbation image block I patch ' and the prior image block C as the guiding condition patch ', the model output is the reconstructed image block I net ; S7: Based on the gold standard image block I in steps S5 and S6 patch , reconstruct image block I net Construct the mean square error loss function L mse and the maximum intensity projection loss function L on three orthogonal planes proj ,constraining the consistency between the model output image and the gold standard image at different observation angles,improving the accuracy of the results; S8: The mean square error loss function L in step S7 is mse With the maximum intensity projection loss L proj Perform weighted calculations and use an optimizer to update parameters of the three-dimensional conditional Poisson generation model; S9: In the model sampling stage, the three-dimensional conditional Poisson generative model trained in step S8 and the prior image block C in step S5 are used. patch , with random noise and C patch As the input of the model, the noise is continuously denoised within a limited step size to obtain the reconstructed image block I net ; will I net The final 3D-DSA reconstructed image is obtained by stitching.

2. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The process of calculating the projection domain image range in step S2 includes constructing the projection matrix, image position point mapping and maximum value calculation. For the point image_index(x,y,z) on the current slice image slice, the rotation angle of the projection image is angle and the voxel spacing is v s , the projection pixel spacing is p s , the distance from the focus to the rotation center is SOD, the distance from the focus to the detector is SID, and the projection matrix M is constructed as Where cos and sin represent the cosine function and sine function respectively, then the scaling coordinate point of the current angle is image_index×v s ×M T , the magnification coeff is SID / (SOD-z), where z represents the z-axis coordinate of the point, The final corresponding projection domain coordinates are point×coeff / p s ; A projection operation is performed on each point of the current slice image, and the maximum value of the projection domain coordinates is calculated, which is used to calculate the range of the projection image that actually participates in training each time, thereby reducing the computing resource consumption during training.

3. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The projection domain image encoder E in step S3 uses U-Net as the backbone network. U-Net is a classic medical image segmentation network. Using U-Net can ensure the model's ability to extract medical image features. The encoder E contains four downsampling modules consisting of convolutional layers, ReLU activation function layers, and maximum pooling layers, and four upsampling modules consisting of convolutional layers and transposed convolutional layers. The input is the processed projection image of N angles, and the output is the corresponding projection image feature f proj The feature conversion process in step 3 includes the construction of the projection matrix, image position point alignment and image interpolation process. The image domain decoder D in step 3 also uses U-Net as the backbone network, and also contains 4 downsampling modules consisting of convolutional layers, ReLU activation function layers and maximum pooling layers, and 4 upsampling modules consisting of convolutional layers and transposed convolutional layers. The input is the converted image feature f slice , the output is the decoded image slice'.

4. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The L1 loss function expression constructed in step S4 is: Where slice is the slice image of each z axis of the gold standard image, and slice' is the image decoded in step S3.

5. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The process of obtaining the prior image C in step S5 is to input N two-dimensional projection images into the model trained in step S4, and then stack each decoded slice' image output on the z-axis; the image cutting process in step 5 is to traverse the coordinates of the corresponding prior image and gold standard image, and cut them with a step size of 64×64×64 and a size of 128×128×128 to obtain the prior image block and the gold standard image block.

6. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The three-dimensional conditional Poisson generative model in step S6 uses 3D-UNet as the backbone network. The training phase of the model includes: (1) Construct random noise of different scales according to data distribution; (2) In each iteration, the gold standard image block I is taken from the gold standard image. patch , add random noise to disturb and get the noisy image block I patch’ ; (3) Noise image block I patch’ and the corresponding prior image block C patch As the input of the 3D conditional Poisson generative model, the output is the reconstructed image block I net .

7. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The maximum intensity projection loss L in step S7 proj The calculation method is to perform maximum intensity projection on the reconstructed image block and the gold standard image block in the x, y, and z planes, and obtain the maximum projection image I of the reconstructed image block in the x-axis direction respectively. net_x , reconstruct the maximum projection image I of the image block in the y-axis direction net_y , reconstruct the maximum projection image I of the image block in the z-axis direction net_z and the maximum projection image I of the gold standard image block in the x-axis direction patch_x , the maximum projection image I of the gold standard image block in the y-axis direction patch_y , the maximum projection image I of the gold standard image block in the z-axis direction patch_z , and then calculate the mean square error loss between them, thereby constraining the consistency of the reconstructed image output by the model and the gold standard image at different observation angles, improving the accuracy of the results, and the maximum intensity projection loss L proj The expression is:

8. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: During model training in step S8, the reconstructed image blocks output by the model and the gold standard image blocks are used to calculate the mean square error loss L MSE and the maximum intensity projection loss L proj To construct the loss function, the comprehensive loss weighted result L total The expression is: THE total =L mse +λL proj Among them, λ is the weighting coefficient, which represents the weight ratio and its value is 0.

1.

9. The sparse view 3D-DSA reconstruction method based on the three-dimensional Poisson generative model according to claim 1, characterized in that: The sampling phase of the three-dimensional conditional Poisson generative model in step S9 includes: (1) Set the maximum scale of the noise, at which the random noise is set as the initial image x 0 ; (2) Set the sampling step size T = 8; (3) For the current step size n, where 0 < n < T, the model input is the current image x n and the prior image patch, and the model output is the image patch x for the next step size n+1 ; then the Euler method is used to solve the ordinary differential equation to correct the result for the next step size; repeat this step until the maximum step size is reached to obtain the final sampling result; (4) The image blocks of each sampling result are spliced ​​according to the spatial position to obtain the 3D-DSA reconstructed image.

Images