A liver tumor radiotherapy dose prediction method based on a diffusion model

Abstract

The application discloses a liver tumor radiotherapy dose prediction method based on a diffusion model, and an embodiment thereof is: 1) data collection; 2) preprocessing; 3) constructing a dose prediction model; 4) constructing a loss function; 5) training the prediction model; and 6) dose prediction. The application introduces a beam field, thereby providing the dose prediction model with radiotherapy beam direction information and dose deposition information in a beam propagation process; in order to efficiently utilize input condition information, in a noise prediction network, a multi-branch encoder and a multi-condition aggregation module are designed to realize feature extraction and aggregation of different input condition information, and an asymmetric fusion module is designed to reduce information loss in the input condition information processing process. Compared with a traditional deep learning method, the application can generate more realistic dose prediction results, is helpful to improve the efficiency of radiotherapy plan design, and formulates a personalized treatment plan, and has great practical application value.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of radiotherapy dose prediction, and particularly relates to a liver tumor radiotherapy dose prediction method based on a diffusion model. BACKGROUND

[0002] Radiotherapy is one of the main means to fight cancer. The design of a radiotherapy plan is the core link of the entire radiotherapy process, and the purpose is to ensure that the irradiation dose received by the tumor target area reaches the required coverage and intensity, and the surrounding normal organs are not irradiated as much as possible. Therefore, dose prediction is a key link in the design of a radiotherapy plan. With the rapid development of deep learning technology, many researchers have begun to explore its application in cancer radiotherapy dose prediction. The existing methods generally use a pixel-level loss function (such as L1 or L2 loss), and the generated dose distribution map is often too smooth, lacking high-frequency details and realism. Although the method based on a generative adversarial network (GAN) can alleviate the problem of the generated result being too smooth to a certain extent, the training process is often not stable enough, and mode collapse is prone to occur, resulting in a lack of diversity in the generated result. In addition, existing researches mainly focus on regions with relatively consistent tumor morphology such as the head and neck and the pelvis, and the difficulty of dose prediction is relatively low. For liver tumors with significant differences in tumor size and location, the prediction complexity of the dose distribution is relatively high, and the existing researches are relatively limited, and further in-depth exploration is needed to improve the prediction accuracy and clinical applicability. SUMMARY

[0003] The present application fully considers the problems existing in the prior art, and aims to provide a liver tumor radiotherapy dose prediction method based on a diffusion model. By optimizing the model structure and parameter configuration, high-precision prediction of the radiotherapy dose of liver tumors is achieved.

[0004] I. Technical principles

[0005] Currently, most dose prediction algorithms take CT images, Planning Target Volume (PTV) contour images, and Organs At Risk (OAR) contour images as model inputs. Due to the significant differences in size and location of liver tumors, it is difficult to effectively train the model to fully capture the complex characteristics of the liver tumor radiotherapy dose distribution, thus making it difficult to achieve accurate dose prediction. To obtain more information to improve the accuracy of dose prediction, the present application constructs a beam field and uses it as a model input, providing the model with beam angle information and dose deposition information during the beam propagation process in liver cancer radiotherapy. To achieve high-precision prediction of liver tumor radiotherapy dose, the present application designs a noise prediction network containing an encoder, an intermediate layer, and a decoder and embeds it into the diffusion model. In the encoder part, feature extraction branches are designed for different input images to achieve targeted feature extraction. To capture the correlation between input images, the present application designs a Multi-condition Aggregation Module (MAM) and uses it in the encoder. Using the Transformer encoder and multi-head attention mechanism, this module can aggregate the features extracted by each branch to extract the dependency between features. In addition, to reduce the information loss caused by image processing, the present application designs an Asymmetric Fusion Module (AFM) and uses it in the decoder, combining the features of the encoder with the features of the decoder to reduce feature loss during input image processing, thereby improving the accuracy of model dose prediction.

[0006] II. According to the above principles, the present application is realized by the following scheme:

[0007] A liver tumor radiotherapy dose prediction method based on a diffusion model, comprising the following steps:

[0008] (1) Data collection

[0009] Collect three-dimensional radiotherapy data of historical liver tumor cases to form a data set, including patient CT images, Planning Target Volume (PTV) contour images, Organs At Risk (OAR) contour images, Intensity-Modulated Radiation Therapy (IMRT) dose maps, and Conformal Radiation Therapy (CRT) dose maps.

[0010] (2) Preprocessing, specifically including the following steps:

[0011] (2-a) Calculate the signed distance map (SDM) based on the PTV contour map and OAR contour map obtained in step (1). The calculation formula is as follows:

[0012]

[0013] Where SDM(u) represents the value at voxel u in the signed distance graph SDM, (i u ,j u ,k u (i) represents the three-dimensional coordinates of voxel u. v ,j v ,k v () represents the three-dimensional coordinates of voxel v, where v represents the surface located in the region of interest (ROI). The voxel points, i.e., the voxel points on the PTV or OAR surface; s i s j s k These represent the voxel spacing in the three different dimensions i, j, and k, respectively; k m It is an integer and k m ∈{-1,0,1},k m =1、

[0014] k m =-1, k m =0 represents voxel u being located inside, outside, and on the surface of the region of interest (ROI), respectively.

[0015] (2-b) Calculate the three-dimensional beam field BF based on the CRT dose map obtained in step (1). The calculation formula is as follows:

[0016]

[0017] Where BF(i,j,k) represents the value of BF at coordinate (i,j,k); CRT(i,j,k) represents the value of the CRT dose map at coordinate (i,j,k), and CRT... min This represents the minimum value in the CRT dose graph. max This represents the maximum value in the CRT dose graph.

[0018] (2-c) Divide the CT image, PTV contour image, OAR contour image, IMRT dose image, CRT dose image, SDM image obtained in step (2-a) and BF image obtained in step (2-b) for each patient in the dataset into two-dimensional images that correspond one-to-one with the patient's original CT slices.

[0019] (3) Constructing a dose prediction model, specifically including the following steps:

[0020] (3-a) The noise adding process gradually adds noise to the IMRT dose map, with the following rules:

[0021]

[0022] where t is the time step, representing different stages in the noise adding process, T represents the maximum time step of noise adding; x t represents the image after adding noise at time step t, x0 represents the IMRT dose map without adding noise; a t represents the noise intensity at time step t, ∈ t is random noise following the standard normal distribution.

[0023] (3-b) The noise removing process gradually removes noise from x T , starting from x T , and finally generates the required dose prediction map, with the following rules:

[0024]

[0025] where x' T = x θ , ∈ t (x′ t , t, y) represents the noise predicted by the noise prediction network from the image x t ' in the noise removing process at time step t; z t is random noise following the standard normal distribution; y represents the conditional information, which is used to guide the model to generate the required dose prediction map, y includes CT images, PTV contour images, OAR contour images, SDM images and BF images.

[0026] (4) Construct the loss function:

[0027] The loss function is constructed as follows:

[0028]

[0029] where, represents the average of the mean square error between the predicted noise ∈ θ (x′ t , t, y) and ∈ t , considering all IMRT dose maps x0 without adding noise, adding noise ∈ t , time step t and conditional information y; MSE(·) represents the mean square error.

[0030] (5) Training the prediction model: input the data obtained in step (2) into the dose prediction model constructed in step (3), train it with the loss function constructed in step (4) as the optimization objective, and update the model parameters using the AdamW optimizer until the loss no longer decreases, to obtain the trained model.

[0031] (6) Dose prediction: use the trained dose prediction model in step (5) to perform dose prediction on the data to be tested, to obtain the dose prediction result.

[0032] The noise prediction network in step (3-b) is composed of an encoder, an intermediate layer, and a decoder.

[0033] The encoder is composed of three branches and a multi-condition aggregation module (MAM).

[0034] Branch 1 includes one 3x3 convolutional layer and five feature extraction modules, and a down-sampling layer is included between every two feature extraction modules, wherein the down-sampling layer includes one 3x3 convolutional layer; the multi-input images composed of the CT image, the SDM image, the PTV contour image, and the OAR contour image are first processed by the 3x3 convolutional layer, and then sequentially pass through the five feature extraction modules to obtain feature maps and X K respectively, and X K is the output of branch 1, wherein the first four feature extraction modules are composed of two residual layers, each residual layer includes two group normalization layers, two Sigmoid linear units, two convolutional layers, an adaptive scaling bias layer, a linear layer, and a residual connection, and the fifth feature extraction module is sequentially stacked by two groups of residual layers-multi-head self-attention layers, each group including one residual layer and one multi-head self-attention layer.

[0035] Branch 3 has the same structure as branch 1, and the BF image is first processed by the 3x3 convolutional layer, and then sequentially passes through the five feature extraction modules to obtain feature maps and X Q respectively, and X Q is the output of branch 3.

[0036] Branch 2 includes one 3x3 convolutional layer and five feature extraction modules, and a pixel-wise addition operation and a down-sampling layer are included between every two feature extraction modules, and the structures of the five feature extraction modules and the down-sampling layer are the same as those of branch 1; the blurred image x t is processed by the 3x3 convolutional layer, and then sequentially passes through the five feature extraction modules, wherein the feature maps generated by the first four feature extraction modules are pixel-wise added with the output feature maps of the corresponding feature extraction modules in branch 1 and branch 3, respectively, to obtain feature maps and Finally, the output feature map X of branch 2 is obtained through the fifth module V .

[0037] The outputs X K , X V , X Q obtained by branches 1, 2 and 3 are input into the multi-condition aggregation module MAM, and the output feature map X E of the encoder is obtained.

[0038] The intermediate layer is sequentially stacked by residual layer 1, multi-head self-attention layer and residual layer 2, wherein the structures of residual layers 1 and 2 are consistent with those of the residual layers in the encoder; the output feature map X E of the encoder is sequentially input into residual layer 1, multi-head self-attention layer and residual layer 2, and the output feature map X M of the intermediate layer is obtained.

[0039] The decoder includes five feature extraction modules and a 3x3 convolution layer, and an up-sampling layer and an asymmetric fusion module AFM are included between every two feature extraction modules, wherein the up-sampling layer includes an interpolation layer and a 3x3 convolution layer, and the structures of the five feature extraction modules are the same as those of the five feature extraction modules of branch 1 in the encoder; the output feature map X M of the intermediate layer is input into the first feature extraction module of the decoder, and the output thereof after the up-sampling layer obtains the feature map D1; D1 and the feature map obtained by branch 2 of the encoder are input into the asymmetric fusion module AFM i , and the output obtained is sequentially input into the second feature extraction module and the second up-sampling layer, and the feature map D2 is obtained; D2 and the feature map obtained by branch 2 of the encoder are input into the asymmetric fusion module AFM2, and the output obtained is sequentially input into the third feature extraction module and the third up-sampling layer, and the feature map D3 is obtained; D3 and the feature map obtained by branch 2 of the encoder are input into the asymmetric fusion module AFM3, and the output obtained is sequentially input into the fourth feature extraction module and the fourth up-sampling layer, and the feature map D4 is obtained; D4 and the feature map obtained by branch 2 of the encoder are input into the asymmetric fusion module AFM4, and the output obtained is sequentially input into the fifth feature extraction module and the 3x3 convolution layer, and the predicted noise ∈ θ of the noise prediction network is obtained. t (x′ t ,t,y).

[0040] The Multi-Condition Aggregation Module (MAM) includes data slicing and flattening operations, positional encoding embedding operations, a Transformer encoder, and linear transformation and reshaping operations. The flattening operation includes a linear layer and a unfolding layer. The positional encoding embedding operation generates randomly initialized positional codes that can be dynamically learned through network training. The output feature maps X from the three branches of the encoder are then used to... Q X K and X V The feature map X' is obtained after slicing and flattening operations and positional encoding embedding operations. Q 、X' K and X' V , will X' Q 、X' K and X' V The inputs are combined with the Transformer encoder to obtain feature map X. Feature map X undergoes linear transformation and reshaping operations to obtain the output feature map X of the MAM. E .

[0041] The asymmetric fusion module (AFM) contains two structurally identical pathways, each including a 3×3 convolutional layer, an instance normalization layer, and a sigmoid unit; the feature map obtained from branch 2 of the encoder is then processed. and the feature map D obtained from the decoder 5-i Input two pathways respectively, where i = 1, 2, 3, 4. Add the feature maps processed by the two pathways pixel by pixel, and then sum them with D. 5-i Asymmetric Fusion Module (AFM) is obtained by splicing together i The output feature map.

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

[0043] (1) This invention introduces the beam field into the model, enabling doctors to interact with the model, which helps to quickly obtain dose prediction maps that meet clinical standards and improve the efficiency and quality of treatment planning.

[0044] (2) The multi-condition aggregation module designed in this invention obtains the data sequence by cutting the data into blocks, thereby using the attention mechanism to capture the relationship between the input condition information and the blurred image at each time step, which helps the noise prediction network generate more accurate prediction results and improves the robustness of the model.

[0045] (3) The asymmetric fusion module designed in this invention integrates features of multi-scale semantic information from the encoder into the decoder, preventing the loss of key information such as the shape, position and boundary of the tumor target area and organs at risk during downsampling, which is conducive to the rapid convergence of the model and improves the accuracy of dose prediction. Attached Figure Description

[0046] Figure 1 A flow chart of a liver tumor radiotherapy dose prediction method based on a diffusion model according to an embodiment of the present application;

[0047] Figure 2 A dose prediction model structure diagram according to an embodiment of the present application;

[0048] Figure 3 A noise prediction network structure diagram according to an embodiment of the present application;

[0049] Figure 4 A multi-condition aggregation module structure diagram according to an embodiment of the present application;

[0050] Figure 5 An asymmetric fusion module structure diagram according to an embodiment of the present application;

[0051] Figure 6 A dose prediction result comparison diagram with other methods according to an embodiment of the present application. DETAILED DESCRIPTION

[0052] The specific embodiments of the present application are described below:

[0053] Embodiment 1

[0054] Figure 1 A flow chart of a liver tumor radiotherapy dose prediction method based on a diffusion model according to an embodiment of the present application is shown, and the specific steps are as follows:

[0055] Step 1, data collection:

[0056] Collect three-dimensional radiotherapy data of historical liver tumor cases to form a data set, including patient CT images, PTV contour images, OAR contour images, IMRT dose images, and CRT dose images.

[0057] Step 2, preprocessing, specifically including the following steps:

[0058] (2-a) Calculate the signed distance map SDM according to the PTV contour image and the OAR contour image obtained in step 1, and the calculation formula is:

[0059]

[0060] Where SDM(u) represents the value at voxel u in the signed distance map SDM, (i u ,j u ,k u ) represents the three-dimensional coordinate value of voxel u, (i v ,j v ,k v) represents the three-dimensional coordinate value of the voxel v, v represents a voxel point located on the surface of the region of interest ROI, i.e., a voxel point on the surface of the PTV or OAR; s i j k m k is an integer and k m ∈{-1, 0, 1}, k m = 1,

[0061] k m = -1, k m = 0 respectively represent that the voxel u is located inside, outside and on the surface of the region of interest ROI.

[0062] (2-b) Calculate the three-dimensional beam field BF according to the CRT dose map obtained in step 1, and the calculation formula is:

[0063]

[0064] Wherein, BF(i,j,k) represents the value of BF at coordinate (i,j,k); CRT(i,j,k) represents the value of CRT dose map at coordinate (i,j,k), CRT min represents the minimum value in the CRT dose map, and CRT max represents the maximum value in the CRT dose map.

[0065] (2-c) Divide the CT map, PTV contour map, OAR contour map, IMRT dose map, CRT dose map, and the SDM map obtained in step (2-a) and the BF map obtained in step (2-b) of each patient in the data set into two-dimensional images corresponding to the original CT slices of the patient one by one, and all images are sampled to be 128*128 in size.

[0066] Step 3, constructing a dose prediction model

[0067] Figure 2 The dose prediction model structure diagram of the embodiment of the application is shown, and the model mainly includes a noise adding process and a noise removing process.

[0068] Noise adding process: gradually add noise to the IMRT dose map, and the rules are as follows:

[0069]

[0070] Wherein t is a time step, which represents different stages in the noise adding process, T represents the maximum time step of noise adding; x t represents the image after noise adding at time step t, x0 represents the IMRT dose image without adding noise; α t ​​​​denotes the noise intensity at time step t, ∈ t is a random noise following standard normal distribution.

[0071] Denoising process: from x T Start to remove noise step by step, finally generate the required dose prediction map, the rules are as follows:

[0072]

[0073] where x' T = x T , ∈ θ (x′ t , t, y) denotes the noise predicted by the noise prediction network from the image x' t in the denoising process at time step t. z t is a random noise following standard normal distribution; y represents conditional information, which is used to guide the model to generate the required dose prediction map, y includes CT image, PTV contour image, OAR contour image, SDM image and BF image.

[0074] Figure 3 The noise prediction network structure diagram of the embodiment is shown, and the network is composed of an encoder, an intermediate layer and a decoder.

[0075] The encoder is composed of three branches and a multi-condition aggregation module MAM;

[0076] Branch 1 includes a 3*3 convolution layer and five feature extraction modules, and a down-sampling layer is arranged between every two feature extraction modules, wherein the down-sampling layer includes a 3*3 convolution layer; the multi-channel input image with a size of 17*128*128 composed of a CT image, an SDM image, a PTV contour image and an OAR contour image is first processed through the 3*3 convolution layer, and then sequentially processed through the five feature extraction modules, to obtain feature maps and X K , X K with a size of 512*8*8 and being the output of branch 1, wherein the first four feature extraction modules are composed of two residual layers, each residual layer includes two normalization layers, two Sigmoid linear units, two 3*3 convolution layers, an adaptive scaling bias layer, a linear layer and a residual connection, wherein the normalization operation adopts group normalization, so that the model can be effectively trained on small batches of data, and the residual connection is combined in the block to avoid the problem of gradient disappearance. The fifth feature extraction module is sequentially stacked by two groups of residual layers-multi-head self-attention layers, each group includes a residual layer and a multi-head self-attention layer with a head number of 4.

[0077] Branch 3 has the same structure as branch 1, and the BF image size is 1x128x128, which is first processed by a 3x3 convolution layer, and then sequentially passes through five feature extraction modules to obtain feature maps and X Q , X Q The size is 512x8x8 and is the output of branch 3.

[0078] Branch 2 includes a 3x3 convolution layer and five feature extraction modules, and a pixel-by-pixel addition operation and a down-sampling layer are included between every two feature extraction modules. The structure of the five feature extraction modules and the down-sampling layer is the same as that of branch 1; the blurred image x t The size is 1x128x128, which is processed by a 3x3 convolution layer, and then sequentially passes through five feature extraction modules, wherein the feature maps generated by the first four feature extraction modules are pixel-by-pixel added to the output feature maps of the corresponding feature extraction modules in branch 1 and branch 3, respectively, to obtain feature maps and The sizes are 128x128x128, 128x64x64, 256x32x32, and 384x16x16, respectively, and finally the output feature map X of branch 2 is obtained through the fifth module V The size is 512x8x8.

[0079] The outputs X K , X V , and X Q obtained by branches 1, 2, and 3, respectively, are input into the multi-condition aggregation module MAM to obtain the output feature map X E of the encoder, which has a size of 512x8x8.

[0080] The intermediate layer is sequentially stacked by residual layer 1, multi-head self-attention layer, and residual layer 2, wherein the structures of residual layers 1 and 2 are consistent with those of the residual layers in the encoder; the output feature map X E of the encoder is sequentially input into residual layer 1, multi-head self-attention layer, and residual layer 2 to obtain the output feature map X M of the intermediate layer, which has a size of 512x8x8.

[0081] The decoder includes five feature extraction modules and a 3x3 convolution layer, and an upsampling layer and an asymmetric fusion module AFM are included between every two feature extraction modules, wherein the upsampling layer includes an interpolation layer and a 3x3 convolution layer, and the structure of the five feature extraction modules is the same as that of the five feature extraction modules of branch 1 in the encoder; the output feature map X MThe first feature extraction module of the input decoder outputs a feature map D1 with a size of 512×16×16 after being upsampled by an upsampling layer. D1 is then compared with the feature map obtained from encoder branch 2. Input Asymmetric Fusion Module (AFM) i The output is then sequentially input into the second feature extraction module and the second upsampling layer to obtain feature map D2, which has a size of 384×32×32. D2 and the feature map obtained from encoder branch 2 are then compared. The inputs are combined with the asymmetric fusion module AFM2, and the outputs are sequentially input into the third feature extraction module and the third upsampling layer to obtain feature map D3, which has a size of 256×64×64. D3 is then combined with the feature map obtained from encoder branch 2. The inputs are fed into the asymmetric fusion module AFM3, and the outputs are sequentially fed into the fourth feature extraction module and the fourth upsampling layer to obtain feature map D4, which has a size of 128×128×128. D4 is then combined with the feature map obtained from encoder branch 2. The inputs are fed into the asymmetric fusion module AFM4, and the output is sequentially passed through the fifth feature extraction module and a 3×3 convolutional layer to obtain the predicted noise ∈ of the noise prediction network. θ (x′ t ,t,y), with a size of 1×128×128.

[0082] Figure 4 The diagram shown illustrates the structure of the multi-condition aggregation module according to an embodiment of the present invention, including data slicing and flattening operations, position encoding embedding operations, a Transformer encoder, and linear transformation and reshaping operations. The flattening operation comprises a linear layer and a unfolding layer. The position encoding in the position encoding embedding operation is randomly initialized and can be dynamically learned through network training. The encoder outputs feature maps X from its three branches. Q X K and X V After the slicing operation, each feature map consists of four non-overlapping blocks. The flattening operation first unfolds each block into a vector of dimension 8192, and then performs a linear transformation to reduce the dimension to 1024. To preserve spatial information, the positional encoding of each block is embedded into the vector to obtain the feature map X'. Q 、X' K and X' V , will X' Q 、X' K and X' V The inputs are combined with the Transformer encoder to obtain feature map X. Feature map X undergoes linear transformation and reshaping operations to obtain the output feature map X of the MAM. E .

[0083] Figure 5 As shown in the structural diagram of the asymmetric fusion module of the embodiment of the application, the asymmetric fusion module comprises two paths with the same structure, each path comprising a 3x3 convolution layer, an instance normalization layer and a Sigmoid unit; the feature map obtained by branch 2 in the encoder and the feature map D 5-i obtained by the decoder are respectively input into the two paths, wherein i = 1, 2, 3, 4, the feature maps processed by the two paths are added pixel by pixel, and D 5-i is spliced to obtain the output feature map of the asymmetric fusion module AFM i .

[0084] Step 4, constructing a loss function

[0085] The loss function is constructed as follows:

[0086]

[0087] wherein, represents the average value of the mean square error between the predicted noise ∈ t (x′ θ , t, y) and ∈ t under the condition of considering all IMRT dose maps x0 without adding noise, adding noise ∈ t , time step t and conditional information y; MSE(·) represents the mean square error.

[0088] Step 5, training a prediction model

[0089] The data obtained in step 2 is input into the dose prediction model constructed in step 3, the loss function constructed in step 4 is used as the optimization objective for training, and the model parameters are updated using the AdamW optimizer until the loss no longer decreases, thereby obtaining the trained model.

[0090] Step 6, dose prediction

[0091] The dose prediction model trained in step 5 is used to perform dose prediction on the data to be tested, thereby obtaining the dose prediction result.

[0092] Example 2

[0093] The method in Example 1 is used to perform dose prediction experiments on a private liver cancer data set. The network implementation is based on PyTorch, and the model is trained on a server equipped with three RTX 2080Ti GPUs. The maximum time step T of adding noise is set to 1000.

[0094] In the present embodiment, the Mean Absolute Error (MAE), HI (Homogeneity Index, HI), CI (Conformity Index, CI), D x x mean max The seven indicators are used to test the U-Net, DeepLabV3+, HD-Net2D, C3D, DoseDiff and MD-Dose six dose prediction networks and the present application on the private liver cancer data set.

[0095] The MAE is defined as the average of the absolute values of the differences between all predicted values and true values, and the calculation formula is:

[0096]

[0097] Where n represents the total number of pixels, represents the true dose value of pixel i. represents the predicted dose value of pixel i. represents the absolute error between the predicted value and the true value.

[0098] D x represents the minimum absorbed dose covering x% of the volume of PTV or OAR. x represents the percentage of the volume of PTV or OAR in which the absorbed dose is greater than or equal to a certain dose threshold x to the total volume of the region. mean represents the average dose absorbed by PTV or OAR. max represents the highest dose absorbed by any voxel in PTV or OAR.

[0099] The results obtained in the experiment are shown in Tables 1 and 2. Table 1 shows the comparison of the MAE indicator, and Table 2 shows the comparison of the HI, CI, D x x mean max It can be found that, compared with other methods, the present application achieves the best result in 16 out of 23 evaluation indicators for PTV and OARs.

[0100] Figure 6 The dose prediction results of the present application and the dose prediction results of other methods are shown in the comparison chart. The dose distribution prediction Figure 1 and the dose distribution prediction Figure 2 The results obtained by the test data are shown in the comparison chart. The difference Figure 1 and the difference Figure 2 represent the dose distribution prediction Figure 1 and the dose distribution prediction​​​​​​Figure 2 The difference with the IMRT manual planning result. The rightmost gray scale bar represents the correspondence between the gray scale in the figure and the numerical value. It can be seen that the prediction error of the method is the smallest, and the beam shape in the prediction result is clear, and the high frequency details are predicted more accurately.

[0101] The above-described embodiments are merely preferred embodiments of the present application and are not intended to limit the scope of the present application. Any changes made to the structure and principles of the present application should be covered within the scope of protection of the present application.

[0102] Table 1

[0103]

[0104] Table 2

[0105]

Claims

1. A method for predicting radiotherapy dose for liver tumors based on a diffusion model, characterized in that... Includes the following steps: (1) Data collection Three-dimensional radiotherapy data of historical liver tumor cases were collected to form a dataset, including patient CT images, contour maps of the tumor target volume PTV, contour maps of organs at risk OAR, dose maps of intensity-modulated radiotherapy (IMRT), and dose maps of conformal radiotherapy (CRT). (2) Preprocessing, which specifically includes the following steps: (2-a) Calculate the signed distance map (SDM) based on the PTV profile and OAR profile obtained in step (1). The calculation formula is as follows: Where SDM(u) represents the value at voxel u in the signed distance graph SDM, (i u ,j u ,k u (i) represents the three-dimensional coordinates of voxel u. v ,j v ,k v () represents the three-dimensional coordinates of voxel v, where v represents the surface located in the region of interest (ROI). The voxel points, i.e., the voxel points on the PTV or OAR surface; s i s j s k These represent the voxel spacing in the three different dimensions i, j, and k, respectively; k m It is an integer and k m ∈{-1,0,1},k m =1, k m =-1, k m =0 respectively indicates that voxel u is located inside, outside, and on the surface of the region of interest (ROI); (2-b) Calculate the three-dimensional beam field BF based on the CRT dose map obtained in step (1). The calculation formula is as follows: Where BF(i,j,k) represents the value of BF at coordinate (i,j,k); CRT(i,j,k) represents the value of the CRT dose map at coordinate (i,j,k), and CRT... min This represents the minimum value in the CRT dose graph. max This represents the maximum value in the CRT dose graph; (2-c) Divide the CT image, PTV contour image, OAR contour image, IMRT dose image, CRT dose image of each patient in the dataset, as well as the SDM image obtained in step (2-a) and the BF image obtained in step (2-b), into two-dimensional images that correspond one-to-one with the patient's original CT slices. (3) Constructing a dose prediction model, specifically including the following steps: (3-a) Noise addition process: Noise is gradually added to the IMRT dose map according to the following rules: Where t is the time step, representing different stages in the noise addition process, and T represents the maximum time step for noise addition; x t This represents the image after adding noise at time step t, where x0 represents the IMRT dose image without added noise; α t Represents the noise intensity at time step t, ∈ t For random noise that follows a standard normal distribution; (3-b) Denoising process, from x T The noise is gradually removed to eventually generate the required dose prediction map, following these rules: Where, x' T =x T ,∈ θ (x t (y, t) represents the noise prediction network's denoising process of image x at time step t. t 'The noise obtained from the prediction;' z t y represents random noise that follows a standard normal distribution; y represents conditional information used to guide the model in generating the required dose prediction map, and y includes CT images, PTV contour images, OAR contour images, SDM images, and BF images. (4) Construct the loss function: The loss function is constructed as follows: in, This indicates that all IMRT dose maps without noise are considered, and those with noise are considered. t Given time step t and conditional information y, the predicted noise ∈ θ (x t ',t,y) and ∈ t The average of the mean squared errors between them; MSE(·) represents the mean squared error; (5) Training the prediction model: Input the data obtained in step (2) into the dose prediction model constructed in step (3), use the loss function constructed in step (4) as the optimization objective for training, and use the AdamW optimizer to update the model parameters until the loss no longer decreases, and obtain the trained model. (6) Dose prediction: The dose prediction model trained in step (5) is used to predict the dose of the data to be tested, and the dose prediction results are obtained.

2. The method for predicting radiotherapy dose for liver tumors based on a diffusion model as described in claim 1, characterized in that... The noise prediction network in step (3-b) consists of an encoder, an intermediate layer, and a decoder; The encoder consists of three branches and a multi-condition aggregation module (MAM); Branch 1 includes a 3×3 convolutional layer and five feature extraction modules. Each pair of feature extraction modules is separated by a downsampling layer, which in turn includes a 3×3 convolutional layer. The multi-channel input image, composed of CT, SDM, PTV, and OAR contour maps, is first processed by the 3×3 convolutional layer, and then sequentially passed through the five feature extraction modules to obtain feature maps. and X K And X K This is the output of branch 1. The first four feature extraction modules consist of two residual layers. Each residual layer includes two groups of normalization layers, two sigmoid linear units, two convolutional layers, one adaptive scaling bias layer, one linear layer, and one residual connection. The fifth feature extraction module consists of two groups of residual layers and multi-head self-attention layers stacked sequentially. Each group includes one residual layer and one multi-head self-attention layer. Branch 3 has the same structure as branch 1. The BF image is first processed by a 3×3 convolutional layer, and then sequentially processed by five feature extraction modules to obtain feature maps. and X Q And X Q This is the output of branch 3; Branch 2 includes a 3×3 convolutional layer and five feature extraction modules. Each pair of feature extraction modules is separated by a pixel-wise addition operation and a downsampling layer. The structure of the five feature extraction modules and the downsampling layer is the same as in Branch 1; (The blurred image x...) t After processing through a 3×3 convolutional layer, the data passes through five feature extraction modules in sequence. The feature maps generated by the first four feature extraction modules are added pixel-by-pixel to the output feature maps of the corresponding feature extraction modules in branches 1 and 3, respectively, to obtain the feature maps. and Finally, the output feature map X of branch 2 is obtained after passing through the fifth module. V ; The outputs X obtained from branches 1, 2, and 3 respectively K X V X Q The multi-condition aggregation module (MAM) is input together to obtain the encoder's output feature map X. E ; The intermediate layer is composed of residual layer 1, a multi-head self-attention layer, and residual layer 2 stacked sequentially, wherein the structures of residual layers 1 and 2 are consistent with those of the residual layers in the encoder; the output feature map X of the encoder is... E The residual layer 1, multi-head self-attention layer, and residual layer 2 are input sequentially to obtain the output feature map X of the intermediate layer. M ; The decoder includes five feature extraction modules and a 3×3 convolutional layer. Between every two feature extraction modules is an upsampling layer and an asymmetric fusion module (AFM). The upsampling layer contains an interpolation layer and a 3×3 convolutional layer. The structure of the five feature extraction modules is the same as that of the five feature extraction modules in branch 1 of the encoder. The output feature map X of the intermediate layer... M The first feature extraction module of the input decoder outputs a feature map D1 after passing through an upsampling layer. D1 is then compared with the feature map obtained in encoder branch 2. Input Asymmetric Fusion Module (AFM) i The output is then sequentially input into the second feature extraction module and the second upsampling layer to obtain feature map D2. D2 is then combined with the feature map obtained from encoder branch 2. The inputs are fed into the asymmetric fusion module AFM2, and the outputs are sequentially fed into the third feature extraction module and the third upsampling layer to obtain feature map D3. D3 is then combined with the feature map obtained from encoder branch 2. The inputs are fed into the asymmetric fusion module AFM3, and the outputs are sequentially fed into the fourth feature extraction module and the fourth upsampling layer to obtain feature map D4. D4 is then combined with the feature map obtained from encoder branch 2. The inputs are fed into the asymmetric fusion module AFM4, and the output is sequentially passed through the fifth feature extraction module and a 3×3 convolutional layer to obtain the predicted noise ∈ of the noise prediction network. θ (x t ',t,y).

3. The method for predicting radiotherapy dose for liver tumors based on a diffusion model as described in claim 2, characterized in that... The Multi-Condition Aggregation Module (MAM) includes data slicing and flattening operations, positional encoding embedding operations, a Transformer encoder, and linear transformation and reshaping operations. The flattening operation comprises a linear layer and a unfolding layer. The positional encoding embedding operation generates randomly initialized positional codes that can be dynamically learned through network training. The output feature maps X from the three branches of the encoder are... Q X K and X V The feature map X' is obtained after slicing and flattening operations and positional encoding embedding operations. Q 、X' K and X' V , will X' Q 、X' K and X' V The inputs are combined with the Transformer encoder to obtain feature map X. Feature map X undergoes linear transformation and reshaping operations to obtain the output feature map X of the MAM. E .

4. The method for predicting radiotherapy dose for liver tumors based on a diffusion model as described in claim 2, characterized in that... The asymmetric fusion module (AFM) contains two structurally identical pathways, each consisting of a 3×3 convolutional layer, an instance normalization layer, and a sigmoid unit; the feature map obtained from branch 2 of the encoder is then processed. and the feature map D obtained from the decoder 5-i Input two pathways respectively, where i = 1, 2, 3, 4. Add the feature maps processed by the two pathways pixel by pixel, and then sum them with D. 5-i Asymmetric Fusion Module (AFM) is obtained by splicing. i The output feature map.

Images