Boron neutron capture therapy Monte Carlo dose super-resolution denoising method

Abstract

The invention discloses a boron neutron capture therapy Monte Carlo dose super-resolution denoising method. A noiseless high-resolution (HR) image is recovered from a noisy low-resolution image (LR). The relationship between the HR image and the LR thereof can be realized through final training output of a CNN-based network, and the effect of MC real-time calculation is achieved through the combination of denoising and super-resolution. According to the scheme, the advantages of high precision of the MC method and accurate and rapid prediction of the CNN are combined, and the problem of time consumption of the traditional MC method is solved. And meanwhile, a physical energy conservation constraint is introduced into network training, so that the problem of non-physical understanding possibly caused by the lack of a physical process of a pure data driven DL method is solved, and metering errors are reduced.

Description

Technical Field

[0001] The invention relates to the technical field of radiotherapy planning systems, in particular to a boron neutron capture therapy Monte Carlo dose super-resolution denoising method. Background Art

[0002] Boron neutron capture therapy (BNCT) has shown great potential in the treatment of head and neck tumors, gliomas, malignant melanomas and other cancers that are poorly treated by traditional radiotherapy methods (such as photon radiotherapy) due to its advantages such as strong targeting and little damage to normal tissues. In recent years, BNCT has developed rapidly in the world, and nearly 20 BNCT projects have been under construction in China by the end of 2024. The therapeutic effect of BNCT directly depends on accurate dose calculation. The Monte Carlo (MC) simulation method is considered to be the "gold standard" for dose calculation. This method can obtain high-precision dose distribution by simulating the transport and various physical processes of neutrons and secondary photons in the patient's body. It is currently the only method that can calculate the dose distribution in three-dimensional inhomogeneous media. However, since MC simulation usually requires billions of particle events to achieve clinically acceptable accuracy (such as uncertainty less than 3%), this process requires hours or even days of computing time, which is difficult to meet clinical scenarios such as online adaptive therapy that require real-time or near-real-time treatment plans.

[0003] Although some studies have used deep learning (DL) technology based on convolutional neural networks (CNN) to suppress Monte Carlo dose noise of photons, no one has applied this technology to BNCT with more complex physical processes. Meanwhile, in image processing, traditional methods focus on single processing of denoising or super-resolution tasks. Pure deep learning methods are prone to produce large dose deviations in key areas (such as skin and tumor interior) in the absence of restrictions on physical processes (such as basic physical laws such as energy conservation), thus affecting the safety and efficacy of treatment. They lack interpretability of the results of clinical dose calculation in BNCT.

[0004] Most of the existing BNCT TPSs use external general Monte Carlo programs such as PHITS and MCNP for dose calculation. Usually, this method divides the calculation area into high-resolution grids and uses a large number of particle event simulations to reduce statistical errors, thereby achieving accurate dose assessment of various organs (including tumors, skin and other key structures). Although this method calculates accurate results, it takes a long time to calculate and has high requirements for hardware resources because it needs to simulate billions of particle events, making it difficult to meet the clinical demand for rapid treatment plan generation.

[0005] Disadvantages of prior art 1: 1. Low computational efficiency: Traditional Monte Carlo programs (such as PHITS and MCNP) are based on CPU architecture, and a single dose calculation takes 6-8 hours (Sumitomo Heavy Industries NeuCure BNCT system in Japan), which cannot meet the real-time requirements of treatment scenarios such as online adaptive radiotherapy (dynamic adjustment of plans during treatment).

[0006] 2. Poor compatibility between hardware and algorithms: General Monte Carlo programs are not optimized for GPU parallel architecture and can generally only run on CPUs. The thread allocation and task batching strategies are extensive, and a large number of computing threads are required to achieve fast calculations, which is not friendly to medical institutions without high-performance computing servers, limiting the clinical use of Monte Carlo software.

[0007] The second existing technology proposes to use a pure deep learning method to achieve BNCT dose prediction. This solution constructs a deep learning model based on 3D U-Net, takes the patient's CT image as input, and directly outputs the three-dimensional dose distribution in the patient's body. By using the nonlinear mapping ability of convolutional neural networks, this method can complete dose calculation in a short time.

[0008] 1.3.2 Disadvantages of the Second Existing Technology 1. Singleness of task: Designed only for denoising tasks, it cannot simultaneously improve spatial resolution and requires additional interpolation steps, resulting in loss of high-frequency details; it cannot further improve the calculation speed 2. Lack of physical constraints: The purely data-driven mean squared error (MSE) loss function ignores physical laws such as energy conservation, and its prediction results often have large deviations in key areas (especially the skin and inside tumors). In addition, the results obtained by pure deep learning predictions are not clinically interpretable enough.

[0009] 3. Poor adaptability of BNCT: The input only contains a single vector diagram, and does not integrate the multi-component doses of multiple secondary particles and patient anatomical information (such as HU value, binary mask of body contour, binary mask of skin area) involved in the BNCT process, making it difficult to model the complex dose superposition effect in BNCT. When encountering special cases that are not fully covered in the training data, unstable prediction results may occur, thus restricting its wide application in clinical practice. Summary of the invention

[0010] The object of the present invention is to provide a boron neutron capture therapy Monte Carlo dose super-resolution denoising method to solve the problems raised in the above background technology.

[0011] To achieve the above object, the present invention provides the following technical solutions: A boron neutron capture therapy Monte Carlo dose super-resolution denoising method comprises the following steps: Step 1: Obtain CT, contouring, and radiotherapy planning data of real patients and use general MC software for dose calculation; Step 2: Divide the acquired data into data sets and perform preliminary processing such as cropping and normalization on the data; Step 3, construct a super-resolution denoising convolutional neural network model; train the traditional super-resolution convolutional neural network model SRCNN and the patented super-resolution denoising convolutional neural network model PDA-SRDNet under the same training conditions to obtain two trained convolutional neural network models; Step 4: Use the test set data to input into two trained convolutional neural networks respectively to obtain the super-resolution denoising results and make multi-dimensional comparisons; Step 5: Perform post-processing operations on the predicted vector map to generate a final vector map; upload the final vector map to the radiotherapy planning system TPS for evaluation by the doctor physicist.

[0012] As a further technical solution of the present invention: In step 1, the general MC program TOPAS is used to set the parameters required for dose calculation, and B-10 (boron dose, D B ), N-14 and H-1 (nitrogen dose and hydrogen dose, D N , D H ) and several absorbed doses and photon doses caused by neutron reactions (D γ Each patient was simulated with 4 doses (D B , D H , D N , D γ ). For each component dose, low-resolution high-noise images and high-resolution low-noise images were simulated as training data and ground truth, respectively.

[0013] As a further technical solution of the present invention: in step 2, the data set is first divided into a training data set and a validation data set. In order to obtain sufficient and diverse training data, each vector map is rotated 90°, 180°, 270°, and flipped up and down, or Gamma mapping is performed, which not only effectively expands the scale of the data set, but also improves the robustness of the network to spatial transformation and noise distribution changes. After data enhancement, we crop, edge fill, and normalize the enhanced data, and uniformly divide it by the maximum value of the low-resolution high-noise vector map.

[0014] As a further technical solution of the present invention: In step 3, the MC super-resolution denoising network selects the improved lightweight convolutional neural network PDA-SRDNet mentioned in this patent, uses four-channel input, and the input includes the patient's CT value three-dimensional matrix, the binary mask of the body contour, the binary mask of the skin area, and the high-noise coarse grid energy deposition matrix calculated by MC; the parameters of the traditional super-resolution network SCRNN and PDA-SRDNet are initialized. The super-resolution denoising network is trained using the binary mask of the body contour, the binary mask of the skin area, and the high-noise coarse grid energy deposition matrix calculated by MC in the training set.

[0015] As a further technical solution of the present invention: in step 4, the comparison of the two convolutional neural networks further includes: performing the same simple processing on the denoising results of the two trained convolutional neural networks obtained in step 3, and after extracting the data, comparing the high-noise low-resolution image obtained from the low-sampling data, the low-noise high-resolution clean image obtained from the high-sampling data, the denoised image of the SRCNN convolutional neural network and the denoised image of the PDA-SRDNet denoising convolutional neural network, and finally outputting the comparison results of the model performance and computational efficiency of the two convolutional neural networks.

[0016] As a further technical solution of the present invention: in step 5, the post-processing operation includes performing an operation of excluding abnormal prediction areas and smoothing edges on the predicted contour. Removing holes, excluding abnormal prediction areas and smoothing edges specifically include the following steps: using a flooding algorithm to remove abnormal hole information within the boundary of the contour; finding the largest connected domain containing the predicted information of the contour to remove other abnormal prediction areas; using the OpenCV built-in library to perform an edge smoothing operation on the contour. The processed denoised image is input into a commercial TPS for evaluation by a physicist.

[0017] Compared with the prior art, the present invention has the following beneficial effects: (1) The image processing acceleration method super-resolution and denoising are effectively combined, and a dual-task joint optimization framework is proposed to accelerate Monte Carlo dose calculation. The high precision of the MC method and the advantages of CNN's accurate and fast prediction are combined to solve the time-consuming problem of the traditional MC method. (2) Introducing strict physical constraints (energy conservation constraints) in network training reduces dosage errors and fundamentally solves the unreliability problem caused by the lack of physical processes in pure DL methods. This is of great significance for the introduction of neural network Monte Carlo denoising into clinical processes.

[0018] (3) Furthermore, a physically constrained super-resolution denoising network consisting of a spatial and channel integrated attention module, a convolutional layer, a pixel shuffle layer, a bilinear, and a residual mechanism is used. While ensuring the denoising performance, a large amount of computational overhead is reduced to reduce the hardware burden, and the patient's dose distribution can be quickly recalculated in a short time, providing instant feedback for adaptive radiotherapy planning, which is suitable for deployment in a clinical real-time computing environment.

[0019] (4) The present invention is not only applicable to dose calculation in BNCT treatment planning, but its basic principle—that is, the statistical denoising method combining MC simulation with super-resolution denoising network—can be extended to other MC radiation transport simulation fields, such as photon therapy planning, proton and heavy ion therapy planning, CT dose assessment, radiation shielding design, and nuclear reactor physics calculation. This method follows a unified statistical and denoising principle and has a wide range of cross-field application potential. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 is an implementation flow chart established in an embodiment of the present invention; Figure 2 is a schematic diagram of network training according to an embodiment of the present invention; Figure 3 Schematic diagram of the PDA-SRDNet network according to an embodiment of the present invention. DETAILED DESCRIPTION

[0021] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0022] See also Figure 1-3 , a Boron Neutron Capture Therapy Monte Carlo dose super-resolution denoising method, which can improve the spatial resolution of the obtained image and remove the inherent noise of the Monte Carlo algorithm. We call this method "Super-resolution Denoising Network in Monte Carlo Denoising Research", and combined with physical constraints, it greatly accelerates the BNCT MC dose calculation while maintaining accuracy, which is of great significance for the clinical application of this method.

[0023] Step 1: Multi-component dose calculation of BNCT based on the general MC software TOPAS: (1) Collect clinical data. Divide the data into data sets: (1.1) Using the general Monte Carlo code TOPAS, the B-10 concentrations in tumor, normal tissue, and skin were set to 60, 20, and 25 ppm, respectively, and there was no B-10 in bone tissue; (1.1.1) Calculate B-10 (boron dose, D) in the same patient model B ), N-14 and H-1 (nitrogen dose and hydrogen dose, D H , D N ) Several absorbed doses caused by reactions with neutrons; (1.1.2) In addition, the primary gamma ray dose from the neutron source and the secondary gamma ray dose produced in the human body by reactions with neutrons and hydrogen make up the photon dose (D γ ); (1.2) Since each absorbed dose component leads to different biological effects, the composite biological effectiveness factor (CBE) and relative biological effectiveness factor (RBE) were used, so that each patient was simulated with a total of 4 doses (D B , D H ,D N , D γ ).

[0024] (1.2.1) According to the latest publication of the International Atomic Energy Agency (IAEA), the “photon isodose” (DisoE) that produces a biological effect equivalent to that of photons can be calculated using formula (1), which is expressed as the weighted sum of the main absorbed dose components:

[0025] in: represents the composite biological effectiveness factor, Indicates the B-10 vector diagram, represents the relative biological effectiveness factor of H, Indicates the H-1 vector diagram; Relative biological effectiveness factor of N, Indicates N-14 vector diagram It shows the relative biological effectiveness factor of Gamma photons; Represents the photon vector diagram.

[0026] (1.3) For the purpose of training, a low-resolution low-sampling noise Monte Carlo simulation was performed on each patient case in the training dataset, with the voxel set to (32×32×24, 9.76×9.76×7.81mm 3 ) Each simulation has 5×10 6For validation purposes, a separate high-resolution, high-sampling “clean” Monte Carlo simulation was performed for each patient case, with a voxel size of (128 × 128 × 96, 2.44 × 2.44 × 1.95 mm 3 ) Each simulation has 1×10 9 neutron events as the true vector diagram.

[0027] (1.4) The reason for predicting high resolution from low resolution is that: although the spatial resolution is reduced, it improves the statistical efficiency. In other words, the result converges faster to the desired particles. This is because a larger sampling intersection surface means that more particle trajectories are scored in each grid, thus reducing the statistical uncertainty. At the same time, from 5×10 6 Vector map prediction 1×10 9 The vector diagram is a process of image denoising.

[0028] (1.5) The goal of using CNN-based methods to solve the denoising and super-resolution problems is to establish a mapping F from the LR image to its HR reconstructed image, which makes the prediction result as similar as possible to the ground truth HR image. Let us denote the LR image as Y. Our goal is to recover the image F(Y) from Y, as shown in Eq. (2). In our proposed method, the goal of PDA-SRDNet is to recover the mapping F of the image F(Y) from Y. The image denoising problem and the SR problem are similar because both problems need to process the high-frequency part while retaining most of the other information. Although the noise in the obtained Monte Carlo vector map is caused by uncertainty, it can still be regarded as a Gaussian noise. Therefore, only one CNN-based deep network is used to solve both problems. This means that there is only one mapping F in our method.

[0029]

[0030] (2) We crop the vector image obtained in step (1) using the body contour as the boundary to eliminate irrelevant areas, such as air, and fill the edges of the valid area (filling each dimension with 0 to the edge of the image).

[0031] (2.1) In order to obtain sufficient and diverse training data, each vector map is rotated 90°, 180°, 270°, flipped upside down, or gamma mapped. This not only effectively expands the size of the dataset, but also improves the network's robustness to spatial transformations and changes in noise distribution.

[0032] (2.2) In order to stabilize the training process, the Monte Carlo dose calculation results are processed. The doses of all patients are normalized by dividing the maximum value of the low-resolution high-noise vector map and linearly normalizing to [0,1]. We calculate it as shown in formula (3):

[0033] (3) Build and train a super-resolution neural network model: Figure 1 As shown in Figure 1, the data preprocessed in steps 1–2 are input into a neural network with an SRCNN structure (or the improved lightweight convolutional network we proposed, such as PDA-SRDNet) to achieve the mapping of low particle number and low resolution dose to high particle number and high resolution dose.

[0034] (3.1) If Figure 2 As shown, the model is constructed using a multi-channel input architecture, and the input includes: a three-dimensional matrix of the patient's HU values; a binary mask of the body contour; a binary mask of the skin area; a low-particle coarse grid energy deposition matrix calculated by MC; (3.2) For the SRCNN network, a simple convolutional layer connection is adopted. SRCNN first uses bicubic interpolation to enlarge the low-resolution image to the target size, and then extracts the features of the image blocks through a three-layer convolutional network. The convolutional layer kernel sizes are 9×9, 1×1, and 5×5 respectively; then the features are nonlinearly mapped and finally a high-resolution image result is reconstructed.

[0035] (3.2.1) The SRCNN network uses MSE as the loss function and uses the Adam learning rate algorithm to train according to the supervised learning algorithm. The formula is as follows:

[0036] in: Represents the network prediction result, represents a low-noise high-resolution vector map, and F represents a mapping from a high-noise low-resolution vector map to a noisy high-resolution vector map.

[0037] (3.3) For the PDA-SRDNet network, Figure 3 As shown in Figure 1, it consists of a spatial and channel integrated attention module, a convolutional layer, a pixel shuffle layer, and a bilinear mechanism. The network contains 6 convolutional layers with dual attention modules, and combines 1 global residual learning (GRL) and 6 local residual learning (LRL).

[0038] (3.3.1) The core component of the network is a lightweight double attention module, which integrates spatial and efficient bidirectional channel attention mechanisms. The spatial attention block captures the pixel-by-pixel relationships within the spatial dimension of feature maps, while the bidirectional channel attention module effectively collects the inter-channel attention information from these feature maps while taking into account both forward and backward dependencies.

[0039] (3.3.2) Residual learning means that the input and output are similar to a large extent, so in our method, the residual is defined as , where most of the values ​​may be zero or very small. In PDA-SRDNet, r is the residual pixel of x and y. Therefore, we can make the whole network training easy and further make the whole network deeper than normal. And both types of this residual, global residual learning and local residual learning, are used in this network. Such a structure can reduce a lot of computational overhead while maintaining accuracy to ensure a balance between denoising accuracy and computational efficiency.

[0040] (3.3.3) The PDA-SRDNet network not only uses the MSE function in formula (4), but also incorporates strict physical constraints into the loss function, introducing physical energy conservation constraints, that is, requiring the energy deposition in each low-resolution unit to be consistent with the sum of the energies of its corresponding N high-resolution units, thereby ensuring that the network prediction satisfies the laws of physics. The loss function used is shown in formula (5). The Adam learning rate algorithm is used for training according to the supervised learning algorithm.

[0041]

[0042] in: Represents the network prediction result, represents a high-noise low-resolution vector map, and N represents the number of noisy high-resolution vector maps corresponding to each voxel of the high-noise low-resolution vector map.

[0043] (4) Super-resolution denoising prediction and dose distribution reconstruction: After training is completed, the high-noise low-resolution dose distribution is input into the trained neural network. The network will quickly output a low-noise fine-grid dose distribution, achieving dose reconstruction from coarse to fine and from high noise to low noise.

[0044] (4.1) and weight the output results of different components according to formula (1) to obtain the final predicted total dose DisoE.

[0045] (4.2) Evaluation indicators include: statistical analysis of the prediction error of each voxel using dose deviation (D_dev); (4.2.1) Use DVH to visually compare dose distributions in critical organs (such as tumors, skin, and other radiation-sensitive tissues); (4.2.2) Calculate the average root mean square error RMSE, GIPR, to quantify the degree of match between the predicted dose and the MC simulation results in terms of space and dose values.

[0046] (4.2.3) GIPR is based on two other dose distribution comparison tools: dose deviation (DD) and distance to agrecment (DTA). Dose deviation uses the dose difference between corresponding points in the reference image and the test image as a measure of similarity. If the dose difference at a point is less than a user-specified value (e.g., 3% of the maximum reference dose), the point is considered to have passed the test. Dose deviation is more suitable for comparing areas with small dose gradients, but is not very reliable in areas with sharp dose changes. The distance to agrecment measures similarity based on the closest distance between points in the reference image and the test image with the same dose. If the distance to agrecment is less than a user-specified value (e.g., 3 mm), the point is considered to have passed the test. In contrast to dose deviation, the distance to agrecment is very efficient in areas with large dose gradients and less efficient in areas with flat doses. The gamma index combines the advantages of both, and its core calculation method is:

[0047] In the formula, and are the coordinates of the midpoints of the reference image and the test image, respectively. For point The gamma index, For point With point The distance between. For point With point The dose difference between . is the threshold of dose deviation, is the threshold of the matching distance. ≤1, then point Passed the gamma test. When calculating the actual gamma pass rate, low-dose areas that are not clinically significant are usually ignored, and the commonly used threshold is 10% of the maximum dose of the reference image.

[0048] (5) Perform post-processing operations on the predicted vector diagram to obtain the final result. Post-processing operations include removing holes, excluding abnormal prediction areas, and edge smoothing to obtain the final dose prediction result, which is uploaded to the radiotherapy planning system TPS for evaluation by the doctor physicist.

[0049] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other specific forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations falling within the meaning and scope of the equivalent elements of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

[0050] In addition, it should be understood that although the present specification is described according to implementation modes, not every implementation mode contains only one independent technical solution. This description of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment may also be appropriately combined to form other implementation modes that can be understood by those skilled in the art.

Claims

1. A method for Monte Carlo dose super-resolution denoising of boron neutron capture therapy, characterized in that: The following steps are involved: Step 1: Obtain CT, contouring, and radiotherapy planning data of real patients and use general MC software for dose calculation; Step 2: Divide the acquired data into data sets and perform trimming and normalization on the data; Step 3, construct a super-resolution denoising convolutional neural network model; train the traditional super-resolution convolutional neural network model SRCNN and the constructed super-resolution denoising convolutional neural network model PDA-SRDNet under the same training conditions to obtain two trained convolutional neural network models; Step 4: Use the test set data to input into two trained convolutional neural networks respectively to obtain the super-resolution denoising results and make multi-dimensional comparisons; Step 5: Perform post-processing operations on the predicted vector map to generate a final vector map; upload the final vector map to the radiotherapy planning system TPS for evaluation by the doctor physicist.

2. The method for Monte Carlo dose super-resolution denoising of boron neutron capture therapy according to claim 1, characterized in that: In the step 1: using the general MC program TOPAS, setting the parameters required for dose calculation, respectively calculating several absorbed doses and photon doses caused by the reaction of B-10, N-14 and H-1 with neutrons on the same patient human body model; simulating a total of 4 partial doses for each patient; and simulating a low-resolution high-noise image and a high-resolution low-noise image for each component dose, respectively, as training data and ground truth.

3. The method for Monte Carlo dose super-resolution denoising of boron neutron capture therapy according to claim 1, characterized in that: In step 2, the data set is first divided into a training data set and a validation data set, each vector map is rotated 90°, 180°, 270° and flipped up and down, or Gamma mapping is performed to achieve data enhancement, the enhanced data is cropped, edge-filled, and data normalized, and uniformly divided by the maximum value of the low-resolution high-noise vector map.

4. The method for Monte Carlo dose super-resolution denoising of boron neutron capture therapy according to claim 1, characterized in that: In step 3, the MC super-resolution denoising network uses an improved lightweight convolutional neural network PDA-SRDNet, and uses a four-channel input, the input including the patient's CT value three-dimensional matrix, the body contour binary mask, the skin area binary mask and the high-noise coarse grid energy deposition matrix calculated by MC; The parameters of super-resolution networks SCRNN and PDA-SRDNet were initialized, and the super-resolution denoising network was trained using the binary mask of the body contour in the training set, the binary mask of the skin area and the high-noise coarse grid energy deposition matrix calculated by MC.

5. The method for Monte Carlo dose super-resolution denoising of boron neutron capture therapy according to claim 1, characterized in that: In the step 4, the comparison of the two convolutional neural networks further includes: performing the same processing on the denoising results of the two trained convolutional neural networks obtained in step 3, extracting the data, comparing the high-noise low-resolution image obtained from the low-sampling data, the low-noise high-resolution clean image obtained from the high-sampling data, the denoised image of the SRCNN convolutional neural network and the denoised image of the PDA-SRDNet denoising convolutional neural network, and finally outputting the comparison results of the model performance and computational efficiency of the two convolutional neural networks.

6. The method for Monte Carlo dose super-resolution denoising of boron neutron capture therapy according to claim 1, characterized in that: In step 5, the post-processing operation includes performing an operation of eliminating abnormal prediction areas and edge smoothing on the predicted contour, specifically including the following steps: using a flooding algorithm to remove abnormal hole information within the boundary of the contour; finding the maximum connected domain containing the predicted information of the contour to remove other abnormal prediction areas; The edge smoothing operation was performed on the contour using the OpenCV built-in library; and the processed denoised image was input into a commercial TPS for evaluation by the physicist.

Images