Four-dimensional cone beam CT reconstruction method based on anatomical prior constraint diffusion model

Abstract

The invention discloses a four-dimensional cone beam CT (Computed Tomography) reconstruction method based on an anatomical prior constraint diffusion model, which mainly solves the problem of poor quality of a reconstructed image in the prior art, and adopts the implementation scheme that the method comprises the following steps: acquiring CT original full-projection data of a batch of patients, and carrying out sparse sampling on the CT original full-projection data; performing data set division and reconstruction on the original full projection data and the sparse projection data to obtain a reference true value and a low-quality CT image containing a large number of artifacts; preliminarily enhancing the CT image containing a large number of artifacts by using the pre-trained ResViT network to obtain a CT image without artifacts; segmenting the CT image without artifacts by using a body segmentation network to obtain anatomical information of each organ; constructing a four-dimensional cone beam CT motion compensation reconstruction model based on anatomical prior constraints, and training and verifying the four-dimensional cone beam CT motion compensation reconstruction model; and inputting the tested data into the trained reconstruction diffusion model to obtain a reconstruction result. According to the method, artifacts of sparse cone beam CT in a specific organ area can be thoroughly eliminated, the reconstructed image is clear in details, and the method can be used for tracking and positioning tumor movement caused by respiration.

Description

Technical Field

[0001] This invention belongs to the field of medical image processing technology, and specifically relates to a four-dimensional cone-beam CT motion compensation reconstruction method, which can be used to track and locate tumor movement caused by respiration. Background Technology

[0002] Cone-beam computed tomography (CBCT) is widely used in image-guided radiotherapy and surgery. CBCT enables precise treatment by monitoring and adjusting patient position in real time, minimizing the impact of patient movement during treatment and improving outcomes. The scan duration is typically about one minute per rotation, covering approximately 10-20 respiratory cycles, which often results in significant motion distortion, particularly in the chest and upper abdomen. These distortions introduce artifacts into the treatment, leading to higher radiation exposure to healthy tissues and a reduced dose to the targeted tumor site.

[0003] Existing four-dimensional cone-beam CT (CBCT) imaging technology can mitigate motion artifacts, encompassing respiratory motion monitoring, projection acquisition, projection classification, reconstruction, and motion artifact correction. It utilizes the patient's respiratory profile to divide multiple projections into different respiratory cycles and reconstructs images from these groups to generate a three-dimensional dynamic sequence of images. However, existing four-dimensional CBCT equipment with respiratory gating devices has long scan times and complex instrument designs, limiting its widespread application. Therefore, while using conventional three-dimensional CBCT scans for four-dimensional reconstruction shows promise, the uneven and sparse angular distribution of each projection group leads to insufficient data, affecting the accuracy of image reconstruction. To address this challenging problem, numerous four-dimensional cone-beam CT reconstruction algorithms have been proposed.

[0004] In their paper "Four-dimensional cone beam ct imaging using a single routine scan via deep learning," published in *IEEE Transactions on Medical Imaging*, Yang et al. proposed a multi-scale adversarial network (MSD-GAN). This method sorts projections based on the diaphragm's position and then reconstructs the sorted data to alleviate severe stripe artifacts in the original image. For the enhanced image, deformation field is used for image registration to estimate the deformation vector field (DVF) between respiratory phases. This method applies the estimated DVF to motion compensation and generates 4D-CBCT images, achieving good results in artifact removal. However, this method only utilizes the local feature extraction capabilities of CNNs, failing to capture global image information. Furthermore, the training process of MSD-GAN is dynamic, requiring continuous interaction between the generator and discriminator, which can lead to instability in the MSD-GAN network.

[0005] Wang et al., in their paper "DuDoTrans: Dual-Domain Transformer for Sparse-View CT Reconstruction" published in *Machine Learning for Medical Image Reconstruction*, proposed a method for post-processing CBCT images while simultaneously restoring the sinusoidal plot. This method uses a deep neural network to interpolate missing data in the sinusoidal plot, thereby recovering and reconstructing the image globally. While this dual-domain approach has achieved widespread success in sparse projection reconstruction tasks and improved image reconstruction quality, the processing of sinusoidal plot data is sensitive; even small changes can significantly impact the reconstructed image and introduce persistent secondary artifacts that are difficult to eliminate. Furthermore, in ultra-sparse view scenarios, interpolation methods cannot accurately restore too much of the original missing data. For example, when interpolating sparse CT data from 40 views to complete CT data from 720 views, the effective information used—the sparse CT data from the 40 views—only accounts for 5.5% of the interpolated data, making learning more difficult and severely impacting reconstruction performance.

[0006] In their paper "DOLCE: A Model-Based Probabilistic Diffusion Framework for Limited-Angle CT Reconstruction," presented at the International Conference on Computer Vision 2023, Liu et al. proposed a method for reconstructing CT images with limited angles based on a diffusion model. This method uses a framework that combines a conditionally trained diffusion model with an explicit physical measurement model to solve the imaging inverse problem. By alternating between data fidelity updates and sampling updates in the conditional diffusion model, high-quality reconstruction of highly ill-posed, angle-limited CT images is achieved. The conditional diffusion model is based on a transformed sine wave. However, because this method reconstructs the entire image using only the CT image and the sine wave, and focuses on specific organs for organ movement caused by respiration, this reconstruction of the entire image may neglect the deformation of specific organs. Summary of the Invention

[0007] The purpose of this invention is to address the shortcomings of the prior art by proposing a diffusion model based on image anatomy prior constraints for sparse projection reconstruction, thereby improving the quality of the reconstructed image.

[0008] The technical approach of this invention is as follows: CT images containing artifacts are directly reconstructed from sparse projection data of cone-beam computed tomography using the Feldkamp-Davis-Kress algorithm. ResViT is used for preliminary data enhancement, and the enhanced images are processed using the embodied segmentation network SwinUNETER to obtain the anatomical structure of organs. The organ anatomical structure is then used as a conditionally driven diffusion model to generate high-quality 4DCBCT images containing information on organ motion caused by respiration.

[0009] Based on the above approach, the implementation steps of this invention include:

[0010] 1. A four-dimensional cone-beam computed tomography (CBCT) reconstruction method based on an anatomical prior constraint diffusion model, characterized in that it includes:

[0011] (1) Obtain a batch of patients’ cone-beam computed tomography raw full projection data from public databases, and sparsely sample them to obtain sparse projection data. Then, preprocess the two types of projection data and divide them into training dataset S1, validation dataset S2 and test dataset S3.

[0012] (2) The Feldkamp-Davis-Kress algorithm was used to reconstruct the original full projection data and sparse projection data to obtain the baseline true value and a low-quality CT image with a large number of artifacts; the ResViT network was pre-trained, and the pre-trained ResViT network was used to perform preliminary enhancement on the CT image with a large number of artifacts to obtain a CT image without artifacts.

[0013] (3) Use the existing embodied segmentation network SwinUNETER to segment CT images without artifacts to obtain anatomical information of each organ;

[0014] (4) Establish the denoising U-Net networks for the forward process, reverse process and reverse process of the diffusion model respectively, and add the anatomical information obtained in step (3) as a condition to the reverse process of the diffusion model to construct a four-dimensional cone-beam CT motion compensation reconstruction model driven by anatomical structure.

[0015] (5) Define the loss function of the diffusion model based on anatomical prior constraints, train the four-dimensional cone-beam CT motion compensation reconstruction diffusion model based on anatomical prior constraints using training dataset S1, and verify it using verification dataset S2 to obtain the trained four-dimensional cone-beam CT motion compensation reconstruction model based on anatomical prior constraints.

[0016] (6) Input the test set S3 into the trained four-dimensional cone-beam motion compensation reconstruction model based on anatomical prior constraints to obtain the reconstruction results of the test set.

[0017] Furthermore, the reconstruction of the original full projection data and sparse projection data using the Feldkamp-Davis-Kress algorithm includes the following implementation:

[0018] 2a) Perform offset correction, gain correction, bad pixel correction and logarithmic transformation on the projected data, and perform cosine weighting on the processed data to compensate for the differences caused by different ray angles and path lengths.

[0019] 2b) For the weighted projection data, perform one-dimensional convolution filtering along the channel direction on each row of detectors;

[0020] 2c) Perform backprojection on the filtered projection, and accumulate the backprojections for each angle of the projection to obtain the reconstruction result;

[0021] 2d) Perform the above operations on the full projection data and sparse projection data respectively to obtain the baseline true value and the low-quality CT image with artifacts.

[0022] Furthermore, step (4) establishes denoised U-Net networks for the forward, reverse, and reverse processes of the diffusion model, respectively, and constructs a four-dimensional cone-beam CT motion compensation reconstruction model based on anatomical prior constraints, the implementation of which includes:

[0023] 4a) By adding noise step by step, the input CT image data is converted into noise data, thus completing the forward propagation process of the diffusion model.

[0024] 4b) Through a series of progressive denoising steps, noise is gradually reduced to reconstruct a high-quality image, thus completing the reverse process of establishing the diffusion model;

[0025] 4c) A U-Net network is used as the denoising network for the inverse process of the diffusion model. It includes an encoder and a decoder, and the feature map of the encoder is passed to the decoder through skip connections to preserve spatial information.

[0026] 4d) The denoising U-Net network, through an additional channel splicing operation between the encoder and decoder, uses the organ anatomical information segmentation map as a spatial constraint to perform deep fusion with the feature map of the current noisy image, ensuring that the reconstructed image can strictly follow the organ morphology and boundary defined by the input anatomical prior while removing noise, thus obtaining a four-dimensional cone-beam CT motion compensation reconstruction model based on anatomical prior constraints.

[0027] Compared with the prior art, the present invention has the following advantages:

[0028] Firstly, this invention uses a whole-body organ segmentation model to extract anatomical information of organs as a conditional driving mechanism.

[0029] The high-quality image reconstruction makes up for the shortcomings of current 4D CBCT image reconstruction that ignores the deformation of specific organs, and can completely eliminate the artifacts of sparse cone-beam CT in specific organ regions.

[0030] Secondly, this invention uses the ResViT network to initially enhance low-quality images and as a conditional input, while simultaneously...

[0031] Global and local conditions provide different scales for image enhancement in the diffusion model, enabling the diffusion model to generate more coherent and natural images.

[0032] Thirdly, this invention ultimately uses a diffusion model as the generative model for generating high-quality 4D CBCT images, and through...

[0033] The overdiffusion model gradually removes noise and refines image details to obtain clear four-dimensional cone-beam CT images, enabling the tracking and localization of tumor movement caused by respiration. Attached Figure Description

[0034] Figure 1 This is a flowchart illustrating the implementation of the present invention;

[0035] Figure 2 This is a schematic diagram of the ResViT network structure used in this invention;

[0036] Figure 3 This is a schematic diagram of the SwinUNETER network structure used to obtain anatomical information in this invention;

[0037] Figure 4 This is a schematic diagram of the reverse process structure of the diffusion model based on anatomical prior constraints in this invention;

[0038] Figure 5 These are comparison images of cone-beam CT images generated by simulation using the present invention and existing methods, respectively. Detailed Implementation

[0039] To enable those skilled in the art to more clearly understand the present invention, the embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the description is only a part of the embodiments of the present invention, and not all of them. Other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort should all fall within the protection scope of the present invention.

[0040] Reference Figure 1 The implementation steps of this example include the following:

[0041] Step 1: Obtain the raw full projection data and perform sparse sampling, preprocessing, and dataset partitioning.

[0042] 1.1) Obtain raw full projection data of cone-beam computed tomography scans of a batch of patients from a public database, and perform sparse sampling on the obtained full projection data to obtain sparse projection data;

[0043] 1.2) The spatial resolution, window level, and window width of the full projection data and sparse projection data are uniformly adjusted, and all data are scaled to a uniform size to obtain the preprocessed full projection dataset S' and sparse projection dataset S.

[0044] 1.3) For each patient in the full projection dataset S' and the sparse projection dataset S, group them according to the respiratory phase to distinguish the respiratory state at different times. That is, derive the respiratory phase information from the data in the sparse projection dataset S according to the diaphragm position and divide each cycle of the respiratory phase into multiple phase segments; then divide all the projection angles contained in the projection into multiple groups of sparse projection data according to different phase segments; divide all the projection angles in the full projection dataset S' according to the same phase segments to obtain the same grouping as the sparse projection dataset S.

[0045] 1.4) Based on the number of patients, the sparse projection dataset S is divided into a 10:1:1 ratio to obtain the training dataset S1, the validation dataset S2, and the test dataset S3. It is ensured that the projection data of the same patient appears in only one of these datasets to prevent inaccurate performance evaluation due to patient data leakage during model training. The full projection dataset S' is also divided in the same way and used as the ground truth for subsequent network training.

[0046] Step 2: Reconstruct the original full projection data and sparse projection data.

[0047] This step uses, but is not limited to, the Feldkamp-Davis-Kress algorithm to reconstruct the sparse projection data and the original full projection data. Its implementation includes the following:

[0048] (2.1) Reconstructing sparse projection data using the Feldkamp-Davis-Kress algorithm:

[0049] 2.1.1) Perform offset correction, gain correction and bad pixel correction on the sparse data in sequence to eliminate the inherent error of the detector itself; then perform logarithmic transformation on the corrected data to convert the X-ray intensity information into line integral projection data that reflects the linear attenuation coefficient of the object.

[0050] 2.1.2) Apply a cosine weighting factor based on cone-beam scanning geometry to the line integral projection data obtained in step 2.1.1). Weighting is applied to compensate for differences in projection values ​​caused by different ray angles and path lengths, resulting in weighted sparse projection data, where:

[0051] The cosine weighting factor Determined by the geometry of the scanning system, its calculation formula is as follows:

[0052] ,

[0053] In the formula The distance from the X-ray source to the center of rotation. and The coordinates of the projection point on the detector plane;

[0054] 2.1.3) For the weighted sparse projection data obtained in step 2.1.2), one-dimensional convolution filtering is performed row by row along the channel direction for each row of data of the detector; the convolution filtering is achieved by convolving the projection data with the ramp filter.

[0055] 2.1.4) The filtered sparse projection data is reconstructed into a 3D image through a backprojection process:

[0056] For each projection angle For each voxel (x, y, z) in three-dimensional space, it is mapped to detector coordinates (u, v) according to the scanning geometry. The corresponding filtered projection value is obtained through interpolation, corrected by a geometric weighting factor, and then accumulated to the voxel. Finally, a CT image containing a large number of artifacts is obtained. Its mathematical expression is as follows:

[0057] ,

[0058] in To reconstruct the 3D image in coordinates The value at that location, The distance from the X-ray to the center of rotation. For the projection angle, In the angle The following is the preprocessed and cosine-weighted projection data. Let v be the impulse response of the ramp filter, and v be the detector row direction coordinate. and For voxel coordinates and scanning angle The intermediate variables calculated based on the system's geometric model are expressed as follows:

[0059] ;

[0060] ;

[0061] 2.1.5) After the processing of steps 2.1.1) to 2.1.4), all sparse projection data are finally converted into low-quality CT images containing a large number of artifacts;

[0062] (2.2) The Feldkamp-Davis-Kress algorithm is used to reconstruct the original full projection data:

[0063] 2.2.1) The original full projection data is sequentially subjected to offset correction, gain correction and bad pixel correction to eliminate the inherent error of the detector itself; then, the corrected data is logarithmically transformed to convert the X-ray intensity information into line integral projection data that reflects the linear attenuation coefficient of the object.

[0064] 2.2.2) The line integral projection data obtained in step 2.2.1) is weighted using a cosine weighting factor based on cone-beam scanning geometry to compensate for differences in projection values ​​caused by different ray angles and path lengths, resulting in weighted full projection data; wherein, the cosine weighting factor... Determined by the geometry of the scanning system, its calculation formula is as follows:

[0065] ,

[0066] in The distance from the X-ray source to the center of rotation. and The coordinates of the projection point on the detector plane;

[0067] 2.2.3) For the weighted full projection data obtained in step 2.2.2), one-dimensional convolution filtering is performed row by row along the channel direction for each row of data of the detector; the convolution filtering is achieved by convolving the projection data with the ramp filter;

[0068] 2.2.4) The filtered full projection data is reconstructed into a 3D image through a backprojection process:

[0069] For each projection angle and each voxel in three-dimensional space Map it to detector coordinates based on the scanning geometry. The corresponding filtered projection value is obtained through interpolation, corrected by a geometric weighting factor, and then accumulated to the voxel to finally obtain the reference true value. Its mathematical expression is as follows:

[0070] ,

[0071] in To reconstruct the 3D image in coordinates The value at that location, The distance from the X-ray to the center of rotation. For the projection angle, In the angle Below is the preprocessed and cosine-weighted projection data. This represents the impulse response of the ramp filter. For the detector row direction coordinates and For voxel coordinates and scanning angle The intermediate variables calculated based on the system geometric model are expressed as follows:

[0072] ,

[0073] ,

[0074] 2.2.5) After the processing described in 2.2.1)-2.2.4), all full projection data are finally converted into high-quality CT images with clear details, which are used as the baseline truth in subsequent steps.

[0075] Step 3: Perform preliminary enhancement on the low-quality image reconstructed by FDK.

[0076] (3.1) Obtain the open-source ResViT network from the publicly available Hugging Face platform:

[0077] Reference Figure 2 The ResViT network structure consists of three cascaded parts, including:

[0078] The first part is the encoder, which uses convolutional layers to transform low-quality CT images containing artifacts into features. The convolutional layers use 3×3 convolutional kernels with a stride of 1 and padding of 1, converting the number of channels of the input image into the required feature dimensions.

[0079] The second part is the information bottleneck module, which consists of an aggregated residual transformer block. The aggregated residual transformer block contains residual connections and a multi-head self-attention structure. The multi-head self-attention layer divides the input feature map into multiple "heads", performs self-attention calculations on each head separately, and then uses the Query, Key and Value matrices to perform self-attention calculations. Finally, the outputs of the multiple heads are concatenated and linearly transformed.

[0080] The third part is a decoder consisting of an upsampling module, which uses a 2×2 convolution kernel with a stride of 2 transposed convolution to progressively restore the spatial resolution of the feature map, doubling the width and height of the feature map, and finally outputting an enhanced CT image of the same size as the original input.

[0081] (3.2) Pre-train the ResViT network to obtain a ResViT network with artifact removal capabilities:

[0082] 3.2.1) Obtain the open-source ResViT network structure from the publicly available Hugging Face platform, and randomly initialize the parameters in the network during the initialization phase;

[0083] 3.2.2) Input low-quality CT images with artifacts into the ResViT network. The network performs nonlinear mapping transformation through multi-layer Transformer modules and residual connections. That is, the network captures long-range dependencies in the image through self-attention mechanism and preserves local detail information through residual structure, thereby gradually learning the complex mapping relationship from image with artifacts to image without artifacts.

[0084] 3.2.3) The output of the ResViT network is compared with the benchmark ground truth image, and the difference between the two is quantified by the designed loss function. The gradient descent optimization algorithm is used to update the parameters in the ResViT network according to the gradient. Through iterative training, the network weights are continuously adjusted so that the model gradually reduces the difference between the output image and the benchmark ground truth, thereby obtaining the ability to remove artifacts.

[0085] 3.2.4) During training, the model performance is evaluated periodically on the validation set, and the optimal model parameters are selected based on the validation loss or image quality index to obtain the ResViT network after pre-training.

[0086] (3.3) The pre-trained ResViT network was used to perform preliminary enhancement on the CT image containing a large number of artifacts:

[0087] The CT image sequence containing a large number of artifacts, obtained by reconstructing sparse projection data using the Feldkamp-Davis-Kress algorithm, is input into a pre-trained ResViT network. Using the prior knowledge learned by the network, the network infers from the input images and directly outputs a CT image sequence without artifacts after preliminary enhancement.

[0088] Step 4: Perform image segmentation on the CT image sequence without artifacts after initial enhancement using the ResViT network to obtain anatomical information of each organ.

[0089] This step employs, but is not limited to, the embodied segmentation network SwinUNETER disclosed by total segmentator to segment artifact-free CT image sequences in order to extract anatomical information of organs from the pre-enhanced CT images. Its implementation includes:

[0090] 4.1) Obtain the pre-trained SwinUNETER network structure and network weights from an open-source website:

[0091] like Figure 3 As shown, the overall structure of the SwinUNETER network follows the U-Net paradigm, consisting of the same number of downsampling and upsampling layers, with skip connections between the feature information extracted from each layer.

[0092] The downsampling layer uses a Swing Transformer block. Each downsampling layer reduces the size of the input feature map to half of its original size and doubles the number of channels in the input feature map.

[0093] The upsampling layer uses transposed convolutional blocks. Each upsampling layer expands the size of the input feature map of the layer to twice its original size and reduces the number of channels of the input feature map of the layer to half its original size.

[0094] After processing the original input image with the same number of upsampling and downsampling layers, the size of the image output by SwinUNETER is the same as the size of the original input image, and the output content is the pixel-level organ segmentation result of the original image;

[0095] 4.2) The CT image sequence without artifacts obtained after the initial enhancement by the ResViT network is used as the input sequence and input into the pre-trained SwinUNETER segmentation network for end-to-end pixel-level semantic segmentation and classification, and the pixels of specific target organs or tissues such as liver and lungs are identified and labeled.

[0096] 4.3) After the above operations, organ anatomical information corresponding to the CT image size is obtained. Each pixel of the input image is assigned a specific category label, which clearly defines the boundaries of the anatomical structures of different organs and tissues. This information is stored in the form of images. These images containing anatomical information serve as prior knowledge for the subsequent construction of motion compensation reconstruction models based on anatomical structures, providing key spatial constraints for the motion compensation reconstruction models.

[0097] Step 5: Construct a four-dimensional cone-beam CT reconstruction model based on an anatomical prior constraint diffusion model.

[0098] This step involves constructing a four-dimensional cone-beam CT reconstruction model based on the fundamental structure of the diffusion model and anatomically constrained diffusion model. This model includes the forward and inverse processes of the diffusion model, as well as the construction of a denoising U-Net network during the inverse process. Its implementation includes the following:

[0099] 5.1) Constructing the forward process:

[0100] The forward process of the diffusion model is a Markov chain whose goal is to progressively degenerate the input data into noise. That is, through a series of steps that progressively add noise, the input CT image data is converted into noisy data. This process is reversible and provides the foundation for the subsequent denoising process, i.e., the inverse process. Its construction method is as follows:

[0101] 5.1.1) Define a standard Gaussian noise distribution in is the standard deviation of the noise, and t is the time step. In this example, the total step size T is set to 1000.

[0102] 5.1.2) Define a noise figure sequence { } is used to control the amount of noise added at each time step; in this example, a linearly increasing sequence is selected.

[0103] 5.1.3) Define the cumulative noise figure This is used to calculate the cumulative noise effect at each time step, the noise figure. The calculation formula is as follows:

[0104] ,

[0105] in It is an iteration of 1-t, representing the index of the time step. It is the noise figure at step i;

[0106] 5.1.4) Add noise gradually, and at each time step, analyze the noise distribution. Medium sampling noise , where I is the identity matrix. The updated noisy image is obtained. :

[0107] ,

[0108] in This is the baseline truth value. The process is a Markov process, where the state at each time step depends only on the state at the previous time step. The image at the current time step t contains noise, and the final image at time step T is obtained. Approaching pure noise;

[0109] 5.2) Constructing the reverse process:

[0110] The reverse process is the opposite of the forward process. It reconstructs the image by gradually removing noise. This reverse process is used to gradually recover the original clear image from noisy data. In this example, the U-Net network is selected as the noise removal network.

[0111] Reference Figure 4The reverse process from the noisy image Initially, for each time step t, the added noise is predicted, thus obtaining further images. :

[0112] ,

[0113] in It is the value of the noise sequence at time step t. It is the noise figure at time step t. It is the noise added during the forward process at the current time step t predicted by the denoising network, and this noise is predicted by the U-Net network in the reverse process;

[0114] 5.3) Construction of the denoising U-Net network during the reverse engineering process:

[0115] Define a U-Net network that includes an encoder and a decoder. The encoder is used to extract features from the input image, and the decoder is used to perform a non-linear mapping on the features extracted by the encoder and output the noise contained in the input image.

[0116] The encoder includes multiple convolutional blocks, each containing a convolutional layer and a pooling layer. The convolutional layer has a kernel size of 3×3 and a stride of 1; the pooling layer has a window size of 2×2 and a stride of 2. The encoder is used to extract information from different levels of the image and fuse it with the information in the corresponding upsampling block.

[0117] The decoder contains multiple upsampling blocks. The number of input blocks to the upsampling blocks is equal to the number of convolutional blocks in the encoder. The upsampling method is transposed convolution with a kernel size of 4×4, a stride of 2, and padding of 1. Each upsampling process doubles the size of the feature map until the output prediction noise of the last layer is the same as the size of the input image of the U-Net network.

[0118] 5.4) U-Net Denoising Network Fusion Prior Analysis:

[0119] Reference Figure 4 The U-Net denoising network uses an additional channel concatenation operation to deeply fuse the organ anatomical information segmentation map as a spatial constraint with the feature map of the current noisy image. This ensures that the reconstructed image, while removing noise, strictly adheres to the organ morphology and boundaries defined by the input anatomical prior. Specifically, the encoder input of the U-Net network is a dual-channel data set. The first channel is the ground truth image after adding noise at the corresponding time step, and the second channel is the anatomical prior information corresponding to the ground truth. The concatenation operation of these two channels deeply fuses the organ anatomical information segmentation map as a spatial constraint with the feature map of the current noisy image.

[0120] Step 6: Train the four-dimensional cone-beam CT motion compensation reconstruction model based on anatomical prior constraints.

[0121] The main training target in this step is the U-Net denoising network used in the inverse process of the diffusion model. The training steps include:

[0122] 6.1) Define the loss function for the reverse process of the diffusion model:

[0123] The output of the denoising U-Net network in the reverse process corresponds to the noise added in the forward process, and its loss function uses the mean squared error loss (MSE):

[0124] ,

[0125] in It represents the number of pixels in the predicted result. It is noise added during the forward process. This is the result predicted by the U-Net network;

[0126] 6.2) Define the hyperparameters for training: Set the learning rate to 10e-4, the batch size to 4, initialize the U-Net network parameters randomly using a Gaussian distribution, and optimize the network parameters using the Adam optimizer;

[0127] 6.3) Training process:

[0128] 6.3.1) Randomly sample artifact-free CT images enhanced by ResViT from the training dataset S1, randomly select a time step, and perform a forward pass to add Gaussian noise of the appropriate level;

[0129] 6.3.2) The noisy image, time-step encoding, and corresponding anatomical segmentation map described in 6.3.1) are input together into the denoising U-Net network, as follows: Figure 4 The schematic diagram of the inverse process of the diffusion model is shown, which obtains the output of the U-Net network, that is, the prediction of the noise added at the current time step;

[0130] 6.3.3) Calculate the mean squared error loss between the network's predicted noise and the actual added noise, backpropagate the loss function value to calculate the gradient of the current U-Net network parameters, and the Adam optimizer uses the gradient of the current U-Net network parameters to update the U-Net network parameters, so that the loss function value continuously decreases.

[0131] 6.4) Periodically validate and save the model:

[0132] 6.4.1) During training, the current model is periodically validated using the validation dataset S2, and its performance on the validation set is calculated. Specifically, the structural similarity (SSIM), peak signal-to-noise ratio (PSNR), and root mean square error (RMSE) of the U-Net-denoised image and the corresponding ground truth image are calculated respectively.

[0133] ,

[0134] ,

[0135] ,

[0136] Where n is the number of observations, i.e. the number of pixels in the calculated image, x is the pixel value of the baseline, and y is the pixel value of the image after denoising by the U-Net network; and y represents the average of the baseline true value x and the pixel value y of the image after denoising by the U-Net network, respectively; and represents the standard deviation of the baseline pixel value x and the standard deviation of the pixel value y of the image after denoising by the U-Net network, respectively. This represents the covariance between the baseline pixel value x and the pixel value y of the image after noise reduction via the U-Net network; and To represent two constants with different values, in order to avoid systematic errors caused by a denominator of zero;

[0137] 6.4.2) Repeat step 6.4.1) When the validation set loss no longer decreases over multiple consecutive training cycles, trigger the early stopping mechanism to terminate the training process and select the model parameters that perform best on the validation set as the final trained four-dimensional cone-beam CT motion compensation reconstruction model to prevent overfitting and ensure the model's generalization ability.

[0138] Step 7: Using the trained four-dimensional cone-beam CT motion compensation reconstruction model, obtain the reconstruction results of the test set.

[0139] The Feldkamp-Davis-Kress algorithm is used to reconstruct CT images with artifacts from the sparse projection data test set S3, and the pre-trained ResViT network saved in step 3 is used to remove the artifacts.

[0140] The SwinUNETER network is used to obtain anatomical information. The image after artifact removal and the anatomical information are then input into the trained four-dimensional cone-beam CT motion compensation reconstruction model. The final reconstruction result is obtained through the denoising network U-Net in the model.

[0141] The effects of this invention can be further illustrated by the following simulation results:

[0142] I. Simulation Conditions

[0143] The simulation platform of this invention is an Ubuntu 22.04 operating system, configured with a 3.50GHz Intel Core i9-10920X CPU and an NVIDIA RTX A5000 GPU, using the PyTorch deep learning framework, developed in Python, version 3.9.

[0144] II. Simulation Content

[0145] Simulation 1: Under the simulation conditions described above, the method of this invention and existing reconstruction algorithms (FDK, Transunet, ResVIT, MSD-GAN, and ResNet) were used to perform high-quality motion-compensated reconstruction on the test set S3 of sparse projected data. The quantized reconstruction results of the three indicators RMSE, PSNR, and SSIM were compared, as shown in Table 1.

[0146] Table 1. Comparison of performance indicators for four-dimensional cone-beam CT reconstruction using different methods.

[0147]

[0148] Among them, RMSE is the root mean square error, PSNR is the peak signal-to-noise ratio, and SSIM is the structural similarity. These three indicators are commonly used image quality assessment metrics. The lower the root mean square error, the higher the peak signal-to-noise ratio and structural similarity, indicating better image quality. As can be seen from Table 1, the method of this invention outperforms other commonly used reconstruction algorithms in terms of overall image reconstruction performance.

[0149] Simulation 2: The data in the test dataset S3 was divided into 10 phases according to the breathing phase. Four phases, phase 1, phase 3, phase 6, and phase 9, were selected. Motion compensation reconstruction was performed using the present invention and existing algorithms such as FDK, Transunet, ResVIT, MSD-GAN, and ResNet. The RMSE, PSNR, and SSIM of each phase were calculated. The results are shown in Table 2.

[0150] RMSE, PSNR, and SSIM indices were calculated to demonstrate the stability of our algorithm across different respiratory phases. Specifically, projection data from five phases were selected, and motion-compensated reconstructions were performed using FDK, Transunet, ResVIT, MSD-GAN, ResNet, and the algorithm of this invention, respectively. The reconstruction results of each algorithm were compared with the baseline values, and the RMSE, PSNR, and SSIM indices were calculated. The final results are shown in Table 2.

[0151] Table 2 Performance metrics of different methods for reconstructed images at different time phases

[0152]

[0153] As can be seen from Table 2, the method of the present invention is stable at different respiratory phases, and its stability at each phase is superior to that of other methods.

[0154] Visualize the parameters in Table 2; the results are as follows: Figure 5 Each column represents a different method, and each row represents the selected respiratory phase: phase 1, phase 3, phase 6, and phase 9, corresponding to the phase rows in Table 2. Meanwhile, Figure 5 Difference maps between the reconstruction results of different methods and the baseline true value were added, located in the four rows of phase 1 difference map, phase 3 difference map, phase 6 difference map, and phase 9 difference map. It can be seen from the difference map that the present invention has a smaller difference from the baseline true value compared with the existing reconstruction algorithms, indicating that the reconstruction effect of the present invention is the best.

[0155] It should be noted that the above steps and the reference numerals in the claims are only for clearly describing the technical solution of the present invention and for ease of understanding, and their order of representation is not limited.

Claims

1. A method for four-dimensional cone-beam CT reconstruction based on diffusion model with anatomical a priori constraints, characterized in that, The method comprises the following steps: (1) obtaining a batch of patient cone-beam computed tomography original full projection data from a public database, sparsely sampling the original full projection data to obtain sparse projection data, and pre-processing the two kinds of projection data to divide the projection data into a training data set S1, a verification data set S2 and a test data set S3; (2) using a Feldkamp-Davis-Kress algorithm to reconstruct the original full projection data and the sparse projection data to obtain a benchmark true value and a CT image with a large number of artifacts of low quality; pre-training a ResViT network, using the pre-trained ResViT network to preliminarily enhance the CT image with a large number of artifacts to obtain a CT image without artifacts; (3) using an existing SwinUNETER segmentation network to segment the CT image without artifacts to obtain anatomical information of each organ; (4) establishing a forward process, an inverse process and a denoising U-Net network of the inverse process of a diffusion model respectively, and adding the anatomical information obtained in step (3) as a condition to the inverse process of the diffusion model to construct a four-dimensional cone-beam CT motion compensation reconstruction model based on an anatomical prior constraint; (5) defining a loss function of the diffusion model, training the four-dimensional cone-beam CT motion compensation reconstruction diffusion model using the training data set S1, and verifying the four-dimensional cone-beam CT motion compensation reconstruction diffusion model through the verification data set S2 to obtain a trained four-dimensional cone-beam CT motion compensation reconstruction model; (6) inputting the test set S3 into the trained reconstruction diffusion model to obtain a reconstruction result of the test set.

2. The method of claim 1, wherein, The pre-processing and data set division of the obtained cone-beam computed tomography original full projection data in the step (1) actually comprise: 1a) a sparse sampling operation: sparsely sampling the full projection data to obtain sparse projection data; 1b) a pre-processing operation: normalizing the full projection data and the sparse projection data, and scaling all the data to a uniform size to obtain a pre-processed full projection data set S' and a sparse projection data set S; 1c) projection data grouping: grouping each patient in the full projection data set S' and the sparse projection data set S according to the respiratory phase to distinguish the respiratory states at different times, 1d) data set division: dividing the sparse projection data set S according to the number of patients in a ratio of 10:1:1 to obtain the training data set S1, the verification data set S2 and the test data set S3, and ensuring that the projection data of the same patient only appears in one of the data sets to prevent inaccurate performance evaluation due to patient data leakage during model training. At the same time, the full projection data set S' is divided in the same way as the benchmark true value for subsequent network training.

3. The method of claim 1, wherein, The reconstruction of the original full projection data and the sparse projection data using the Feldkamp-Davis-Kress algorithm in the step (2) comprises the following steps: 2a) performing offset correction, gain correction, bad point correction and logarithmic transformation on the projection data to prepare for subsequent processing; 2b) performing cosine weighting on the data processed in step 2a) to compensate for the differences caused by different ray angles and path lengths; 2c) One-dimensional convolution filtering is performed on the weighted projection data along the channel direction on each row of the detector; 2d) The filtered projection is back-projected, and the back-projection is accumulated for each angle of the projection to obtain the reconstruction result; 2e) The above operations are performed on the full projection data and the sparse projection data respectively to obtain the benchmark true value and the CT image with low quality and artifacts.

4. The method of claim 1, wherein, The pre-training of the ResViT network in (2) includes: 2f) Obtain the open-source ResViT network structure from a public website, input the low-quality image with artifacts into the ResViT network, perform nonlinear mapping transformation through the Transformer layer network structure in the ResViT network, compare the output result of the ResViT network with the benchmark true value, obtain the difference between the two, and perform back propagation to continuously correct the parameters in the ResViT network, so that the ResViT network obtains the ability to remove artifacts, and a pre-trained ResViT network is obtained.

5. The method of claim 1, wherein, In (2), the pre-trained ResViT network is used to preliminarily enhance the CT image containing a large number of artifacts. The CT image sequence containing a large number of artifacts obtained by reconstructing the sparse projection data using the Feldkamp-Davis-Kress algorithm is input into the pre-trained ResViT network. The network uses the prior knowledge it has learned to infer the input image and directly outputs a CT image sequence without artifacts.

6. The method of claim 1, wherein, In (3), the existing SwinUNETER embodiment segmentation network is used to segment the CT image without artifacts, which includes: 3a) Obtain the pre-trained SwinUNETER network structure and network weight from an open-source website; 3b) Input the CT image sequence without artifacts obtained after preliminary enhancement by the ResViT network into the pre-trained SwinUNETER segmentation network to perform pixel-level classification on the organ information in the CT image; 3c) The SwinUNETER network outputs segmentation maps corresponding to the size of the input image, and extracts the accurate spatial profile and position information of the target organ from the segmentation maps as the prior constraint for constructing the motion compensation reconstruction model based on anatomical prior constraint.

7. The method of claim 1, wherein, In (4), the forward process, reverse process, and denoising U-Net network of the diffusion model are established, and a four-dimensional cone beam CT motion compensation reconstruction model based on anatomical prior constraint is constructed, which includes: 4a) Construct the forward process of the diffusion model: The forward process of the diffusion model is a Markov chain, and the goal is to gradually degrade the input data into noise, that is, through a series of steps of gradually adding noise, the input CT image data is converted into noise data; 4b) Construct the reverse process of the diffusion model: The goal of the reverse process of the diffusion model is to gradually recover the original clear image from the noisy data, that is, through a series of steps of gradually denoising, the noise is gradually reduced to reconstruct a high-quality image; 4c) Construct the denoising U-Net network of the reverse process of the diffusion model: The U-Net network is used as a denoising network of the inverse process of the diffusion model, which comprises an encoder and a decoder, and the feature maps of the encoder are transmitted to the decoder through a skip connection to retain spatial information; The encoder comprises a plurality of convolution blocks, each of which comprises a convolution layer and a pooling layer, for extracting information of the image; The decoder comprises a plurality of upsampling blocks, the number of the upsampling blocks is equal to the number of the convolution blocks in the encoder, and the upsampling mode is transposed convolution; 4d) Network structure fusion of anatomical prior: The denoising U-Net network fuses the organ anatomical information segmentation atlas as a spatial constraint with the feature maps of the current noisy image through an additional channel splicing operation between the encoder and the decoder, so as to ensure that the reconstructed image can strictly follow the organ morphology and boundary defined by the input anatomical prior while removing noise, and obtain a four-dimensional cone beam CT motion compensation reconstruction model based on the anatomical prior constraint. The loss function of the diffusion model defined in (5) is used to train the four-dimensional cone beam CT motion compensation reconstruction diffusion model by using the training data set S1, and the model is verified by using the verification data set S2, and the implementation includes:

8. The method of claim 1, wherein, 5a) Defining the loss function of the inverse process of the diffusion model: The output of the denoising U-Net network of the inverse process is the noise added in the forward process, and the loss function uses the mean square error loss MSE: 5b) Randomly sampling a clean CT image, i.e., a reference true value, from the training data set S1, randomly selecting a time step, adding a corresponding level of Gaussian noise by executing the forward process, and inputting the noise image, time step code and corresponding anatomical segmentation atlas into the denoising U-Net network; calculating the mean square error loss of the network predicted noise and the real added noise, and updating the network parameters by using the gradient descent method iteratively; ； wherein is the number of predicted result pixels, is the noise added in the forward process, is the result predicted by the U-Net network; 5c) During the training process, the verification data set S2 is periodically used to verify the current model, and the loss on the verification set is calculated; when the loss on the verification set no longer decreases in continuous multiple training periods, the early stopping mechanism is triggered, the training process is terminated, and the model parameters with the best performance on the verification set are selected as the final trained four-dimensional cone beam CT motion compensation reconstruction model, so as to prevent overfitting and ensure the generalization ability of the model. ​

Images