CN107551411B - Proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty - Google Patents
Proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty Download PDFInfo
- Publication number
- CN107551411B CN107551411B CN201710707207.6A CN201710707207A CN107551411B CN 107551411 B CN107551411 B CN 107551411B CN 201710707207 A CN201710707207 A CN 201710707207A CN 107551411 B CN107551411 B CN 107551411B
- Authority
- CN
- China
- Prior art keywords
- range
- heavy ion
- dose
- proton heavy
- radiotherapy
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Radiation-Therapy Devices (AREA)
Abstract
The invention discloses a robust optimization method of proton heavy ion intensity modulated radiation therapy aiming at range uncertainty. The method not only can accurately establish the range uncertainty model of the robust optimization problem, but also has the characteristics of high solving speed, stable algorithm, strong robustness and the like, so that a user can reduce the influence of range errors on a radiotherapy plan through the robust optimal solution.
Description
Technical Field
The invention belongs to the field of optimization methods of radiotherapy, and relates to a proton heavy ion modulated intensity radiotherapy robust optimization method aiming at range uncertainty.
Background
Proton heavy ion radiotherapy technology is one of the most advanced radiotherapy means, can effectively protect normal tissues (higher treatment gain ratio) while accurately killing tumor cells, and is a 'sharp tool' for tumor treatment. When the proton heavy ion radiotherapy technology reaches a tumor focus, a large amount of energy is released instantly to kill tumor cells, and then the energy is attenuated rapidly, so that the dosage is called as a Bragg peak in the dosage science, and the pre-peak dosage is about 20 percent of the peak value. Proton heavy ion modulated radiotherapy obtains the expected dose targets of a target area and a endangered organ by optimizing the energy and the intensity of a proton heavy ion scanning beam, is the main mode of proton heavy ion radiotherapy, and has great influence on the precision by range, respiratory motion and the like. Proton heavy ion radiotherapy is like 'three-dimensional directional blasting', and has the characteristics of effectively killing tumors without generating obvious toxic and side effects. Therefore, compared with photon radiotherapy, proton heavy ion radiotherapy has more advantages in treating children tumors and tumors (chordoma, chondrosarcoma and the like) close to important organs, and can obviously reduce the occurrence probability of secondary primary tumors.
The treatment precision of proton heavy ion modulated radiotherapy is high, so error management is very important. Errors (namely uncertainties) can exist in multiple links in proton heavy ion modulated radiotherapy, and the errors cannot be eliminated and can only be reduced. If the proton heavy ion beam bragg peak depth is below or above the expected depth, organs at risk near the tumor may be exposed to high doses, causing complications of radiotherapy, and even life-threatening. Unlike photon radiotherapy, proton heavy ion modulated radiotherapy has difficulty in reducing the effects of uncertainty by way of tumor target volume (GTV) flaring to form a Planned Target Volume (PTV). Although proton heavy ion modulated radiation therapy can produce a better dose distribution than photon radiation therapy, it is difficult to further improve the efficacy if the uncertainty is not effectively controlled.
Research generally considers that one of the main factors influencing the accuracy of proton intensity modulated radiotherapy is range uncertainty. Expert scholars propose two main sources of uncertainty of proton intensity modulated radiotherapy: the range uncertainty caused by CT artifact, field unit conversion, anatomical structure change and the like; the second is patient modeling uncertainty due to patient position deviation or tissue density non-uniformity. There are experts that believe that proton range uncertainty comes from a variety of sources, mainly due to patient geometric anatomical variations, and in addition treatment physics and calculations also contribute to range uncertainty.
although range uncertainty of proton heavy ion intensity modulated radiotherapy cannot be avoided, negative effects of proton heavy ion intensity modulated radiotherapy can be reduced by a robust optimization range uncertainty model. At present, the uncertainty model of proton heavy ion intensity modulated radiotherapy range mainly comprises two types: a maximum minimum model and a worst dose distribution model. The maximum minimum model formula is as follows:
minx{maxA'∈R{f(g(x,A'))},x≥0} (5)
Wherein x ∈ RnIntensity vector representing scanning point, A' is belonged to R and represents uncertainty solution set, R is { A }1,A2,…,AnIs the dose contribution matrix uncertainty solution. The function g represents the dose for all voxel points on the organ, g (x, a ') ═ a' x. f is a penalty function representing the deviation between the optimized dose and the desired dose. The range of the proton heavy ion beam is considered to be a random variable. Therefore, the irradiation dose and the objective function value are also random variables. A general randomness model is built by introducing uncertainty solution probability p (a'), such as:
minx{max{∑A'∈Rp(A')·f(g(x,A'))},x≥0} (6)
Wherein p (A'). epsilon.R|R|，∑A'∈Rp (a') is defined empirically as 1. As can be seen from equation (6), the maximum-minimum modelThe worst solution is found out from the range uncertainty solution set and then is minimized through optimization.
The worst dose distribution model requires statistics of the minimum tumor target dose and the maximum organ-at-risk dose in the range uncertainty interval. Although the worst dose distribution model does not physically hold, a conservative lower bound is provided. The model formula is as follows:
wherein A is0A solution set representing the desired value is presented,And λA'Weight factors, g, representing the expectation function and worst dose penalty function, respectivelynormalRepresents the dose of the voxel under normal conditions, gworst-caseRepresents the dose at the voxel point in the worst case and f is a penalty function.
The main problem with both types of models described above is that it is not reasonable to obtain a robust optimal solution based only on the worst "bad" case of the current solution within the perturbation interval.
Particle Swarm Optimization (PSO) is a heuristic proposed in 1995 to enable a swarm to achieve foraging purposes through collective cooperation and competition of a bird swarm. Definition f (x) is in D dimension Euclidean space EDIs a real function of. x is the number ofi＝(xi1,xi2,…,xiD)Tindicating the current position of the ith particle. p is a radical ofi＝(pi1,pi2,…,piD)TIndicating the optimal position, p, that the ith particle has undergoneg＝(pg1,pg2,…,pgD)TRepresenting the best position of the population. v. ofi＝(vi1,vi2,…,viD)TIndicates the speed, fitness, of the ith particlei＝f(xi) The fitness value of the ith particle is expressed. The particles generate a new generation of individuals in a speed and position updating mode, see formulas (8) and (9), and the whole population is continuously iterated to achieve the purpose of optimizing.in general, the change in position of the particle in the D (1. ltoreq. D. ltoreq. D) th dimension is limited to a regionWithin, the speed variation is limited to regionsAnd (4) the following steps.
vid(t+1)＝ω·vid(t)+c1·r1·(pid(t)-xid(t))+c2·r2·(pgd(t)-xid(t)) (8)
xid(t+1)＝xid(t)+vid(t+1) (9)
where i is 1,2, …, M, D is 1,2, …, D, M is the population number, ω is the inertial weight, c is the inertial weight1As a cognitive parameter, c2Is a social parameter, r1And r2Are respectively two mutually independent in [0,1 ]]A random number in between.
The particle swarm optimization has good performance in solving the robust optimization problem, but the particle swarm optimization problem is easy to fall into a local extreme value.
Disclosure of Invention
In order to solve the defects existing in the prior art, the invention provides a proton heavy ion modulated intensity radiotherapy robust optimization method aiming at range uncertainty, an uncertainty model of a robust optimization problem can be accurately established, the model is solved by adopting a self-adaptive particle swarm robust optimization method based on a filter algorithm, the solving speed is high, the algorithm is stable, the robustness is strong, and a user can obtain a robust optimal solution.
The invention adopts the following technical scheme for solving the technical problems:
The proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty comprises the following specific steps:
Step 1: constructing a range uncertainty mathematical model based on the combination of optimality and robustness; the model embodies the optimality and robustness of radiotherapy plans caused by range uncertainty, and specifically comprises the following steps:
In the above-mentioned formula,The objective function values at different ranges are represented, and optimality is reflected.The robustness is embodied, namely the deviation of the target function value under different ranges. w represents the intensity of the pencil beam of the proton heavy ion beam, z represents the organ sampling point at rangeThe actual dose to be administered.represents a proton heavy ion beam set, B (i) represents a pencil beam set of beam i (i ∈ B). Range of jth pencil beam of ith beamIt is shown that,Indicating the normal range of the gas flow in the gas flow path,Indicating a maximum above/below normal range, the range of the jth pencil beam of the ith beam is indicated as For an arbitrary beam i ∈ B, use is made ofRepresenting the set of possible ranges after uncertainty dispersion of pencil beam j ∈ b (i), and R represents the number of elements in the range set. Pencil beams j andEach element in the array is in one-to-one correspondence, so that the array can be at any rangeby indexingand (6) unique confirmation. Use ofIndex element representing rangeThe vector of (2). For convenience of definitionIs at normal rangeIs used to determine the index of (1). By usingrepresenting the Target region (Target) and the Organ At Risk (OAR), respectively, and v represents the v th voxel.AndRepresenting the weights of the target and the organs at risk, respectively. By wjIndicating the intensity, theta, of the jth pencil beamvRepresenting the desired dose value for the ν th voxel.
Step 2: and optimizing by adopting a self-adaptive particle swarm robustness optimization method based on a filtering algorithm to obtain a robust optimal solution.
Step 2.1: the solution violating the constraint is first "filtered" using a filtering algorithm. The solution in equation (1) represents the weight of pencil beams, the weight is a non-negative number, and if the weight obtained in the optimization process is negative, the processing is performed in a filtering mode of setting zero.
Step 2.2: in order to avoid the particle swarm optimization from falling into local optimization, an effective dynamic variation strategy is adopted. Introducing a diversity function D and a fitness variance function sigma2。
f＝max{1,max{|f(i)-favg|}} 1≤i≤M
wherein M is the population number and L is the search radius. f (i) is the fitness function value of the ith particle, favgAnd the average value of the fitness function values of the current particle swarm is obtained.
the variation function of particle velocity is shown below:
Wherein v isiis the velocity of the ith particle, t is the number of iterations, GmaxFor the total number of iterations, θ is a random number that follows a standard normal distribution N (0, 1).
Step 2.3: optimizing formula (1) by adopting a particle swarm algorithm, filtering solution violating the constraint by adopting the method in step 2.1 in the optimization process, and calculating the particles by adopting formula (2) and formula (3) to obtain D and sigma2If D and σ2And if the threshold value exceeds the preset threshold value, performing mutation operation by adopting a formula (4). And after the stopping criterion is met, finally obtaining a robust optimal plan of the proton heavy ion intensity modulated radiotherapy aiming at the range uncertainty.
the method further comprises the following sequence operations before the step 1:
1) importing radiotherapy plan information
Radiotherapy planning information comprises CT images of a patient, beam parameters, accelerator hardware parameters, organ constraint parameters and range uncertainty parameters; a standard DICOM RT format is used as an input and output interface, a radiotherapy plan file of a commercial plan system is imported, after a radiotherapy plan is finished, proton heavy ion scanning point weight and dose information are exported to the DICOM file, and the supported standard files mainly comprise: CT, DICOM RT-Structure, DICOM RT-Plan, DICOM RT-Dose;
2) Calculating a dose contribution matrix
calculating the dose contribution of proton heavy ion scanning points with different energies and different positions to all organ sampling points, adopting a fast analysis method-pencil beam dose calculation algorithm, and performing convolution superposition on the flux weight of the proton heavy ion at the position to be calculated and IDD data corresponding to the water equivalent depth, wherein the formula can be expressed as follows:
Dose(r，z)＝ΣkΦk(r，z)·IDD(z) (10)
Wherein z represents the depth position along the beam, r represents the transverse position, k represents the index of the proton heavy ion pencil-shaped beam spot, Dose represents the total Dose in a voxel corresponding to a certain position, phi is the flux weight of the proton heavy ion, and IDD is the integral depth Dose data of the proton heavy ion pencil-shaped beam in the depth direction in water which is calculated in advance through actual measurement or a Monte Carlo method.
The sequence of the step 2) comprises the following specific steps:
Firstly, distributing points according to the full width at half maximum of a proton heavy ion beam spot in the beam direction; dividing each beam spot into 19 sub-beams; then tracking each sub-beam by adopting an improved Siddon algorithm, calculating the water equivalent depth corresponding to each voxel on a ray path, and looking up a table to obtain a corresponding IDD value; according to the Fermi-Eyges theory, a two-dimensional calibration mode is adopted to calculate the Gaussian distribution standard deviation sigma of the Torontal scattering corresponding to the voxel passed by the rayMCS(ii) a Then, the influence of nuclear scattering on beam flux transverse scattering is calculated according to a theoretical formula provided by Soukup and the like to obtain the weight omega of the nuclear scatteringNSAnd the standard deviation sigma of its Gaussian distributionNS(ii) a Considering the truncation error δ of the gaussian distribution, the total flux weight formula is obtained as:
Φ(x,y,z)＝(1-ωNS)·[G(x,y,σMCS)+δMCS]+ωNS·[G(x,y,σNS)+δNS] (11)
According to beam monitoring system sensitivity factor SmAnd superposing all the results according to weight to obtain the final dose (x, y, z):
dose(x,y,z)＝∑e∑jSm·ωe,j·De,j(x,y,z) (12)
Where x, y, z represent position coordinates, e represents the energy layer index, j represents the pencil beam spot index, SmIs a sensitive factor, omega, of the monitoring system related to the beam energye,jis the weight of the jth beam of energy layer e, De,j(x, y, z) is the dose of energy layer e jth beam vs. position (x, y, z).
Compared with the prior art, the invention has the beneficial effects that:
1. The invention constructs a novel range uncertainty model of proton heavy ion intensity modulated radiotherapy, and overcomes the defect that the prior model is optimized according to the worst condition of the current solution in a disturbance interval; in addition, according to the characteristics of a range uncertainty model of proton heavy ion intensity modulated radiotherapy, a self-adaptive particle swarm robustness optimization method based on a filter algorithm is provided, the robustness optimization effect is improved, and the treatment effect of a proton heavy ion radiotherapy plan is improved.
2. The novel range uncertainty model constructed by the method is more reasonable. The maximin model and the worst dose distribution model researched by the predecessor can reduce the negative influence of range uncertainty of proton intensity modulated radiotherapy to a certain extent, but only the worst 'difference' condition of the current solution is considered, and the distribution of the solution in an uncertainty interval is not considered, so that the solved solution is over-conservative and the robust optimal solution is difficult to obtain; the method establishes a range uncertainty model which simultaneously embodies optimality and robustness, and improves the reasonability and the practicability of the model.
3. The method for solving the range uncertainty model by adopting the adaptive particle swarm robustness optimization method based on the filtering algorithm is more effective and faster; although early domestic and foreign research also involved robust optimization, it was basically based on improvements of conventional optimization methods such as L-BFGS-B code and CPLEX toolkit; the conventional optimization method has strong universality, but the uncertainty problem solving capability is general, so that the optimization time is long (the longest time can be tens of hours), the optimization effect is to be improved, and the clinical application of the method is restricted; therefore, the invention has extremely important medical application value and overcomes the unsolved technical problem.
Drawings
in order to facilitate understanding for those skilled in the art, the present invention will be further described with reference to the accompanying drawings.
FIG. 1 is an overall flow chart of the present invention;
Fig. 2 is a reference geometric anatomical model of the present invention.
Detailed Description
According to the embodiment, proton heavy ion intensity modulated radiotherapy plan information and preparation calculation data are firstly introduced, then a range uncertainty model based on combination of optimality and robustness is constructed according to a robustness measurement method of mean and variance, and finally a robust optimal solution is obtained by adopting adaptive particle swarm robust optimization based on a filtering algorithm. The method can accurately establish the uncertainty model of the robust optimization problem, has high solving speed, stable algorithm and strong robustness, and can ensure that a user can obtain the robust optimal solution.
as shown in fig. 1, the present invention provides a robust optimization method for proton heavy ion modulated radiotherapy with uncertainty of range, which comprises the following steps:
(1) data preparation
1) importing radiotherapy plan information
Radiotherapy planning information mainly comprises patient CT images, beam parameters, accelerator hardware parameters, organ constraint parameters, range uncertainty parameters and the like. The project adopts a standard DICOM RT format as an input and output interface, can import a radiotherapy plan file of a commercial plan system, and export the weight and the dose information of proton heavy ion scanning points to the DICOM file after completing the radiotherapy plan, wherein the supported standard files mainly comprise: CT, DICOM RT-Structure, DICOM RT-Plan, DICOM RT-Dose.
2) Calculating a dose contribution matrix
The dose contribution of proton heavy ion scanning points with different energies and different positions to all organ sampling points needs to be calculated, the calculated amount is large, and a relatively mature rapid analysis method, namely a pencil beam dose calculation algorithm, is adopted. The core idea of the method is a convolution superposition method, and specifically, the flux weight of the proton heavy ion at the position to be calculated and the IDD data corresponding to the water equivalent depth are subjected to convolution superposition. Can be expressed by the formula:
Dose(r,z)＝∑kΦk(r,z)·IDD(z) (10)
Wherein z represents the depth position along the beam, r represents the transverse position, k represents the index of the proton heavy ion pencil-shaped beam spot, Dose represents the total Dose in a voxel corresponding to a certain position, phi is the flux weight of the proton heavy ion, and IDD is the integral depth Dose data of the proton heavy ion pencil-shaped beam in the depth direction in water which is calculated in advance through actual measurement or a Monte Carlo method.
Firstly, distributing points according to the full width at half maximum of a proton heavy ion beam spot in a beam direction view (BEV) direction; dividing each beam spot into 19 sub-beams; then tracking each sub-beam by adopting an improved Siddon algorithm, calculating the water equivalent depth corresponding to each voxel on a ray path, and looking up a table to obtain a corresponding IDD value; according to the Fermi-Eyges theory, a two-dimensional calibration mode is adopted to calculate the Gaussian distribution standard deviation sigma of the Torontal scattering corresponding to the voxel passed by the rayMCS(ii) a Then, the influence of nuclear scattering on beam flux transverse scattering is calculated according to a theoretical formula provided by Soukup and the like to obtain the weight omega of the nuclear scatteringNSAnd the standard deviation sigma of its Gaussian distributionNS(ii) a Considering the truncation error δ of the gaussian distribution, the total flux weight formula is obtained as:
Φ(x,y,z)＝(1-ωNS)·[G(x,y,σMCS)+δMCS]+ωNS·[G(x,y,σNS)+δNS] (11)
According to beam monitoring system sensitivity factor Smand superposing all the results according to weight to obtain the final dose (x, y, z):
dose(x,y,z)＝∑e∑jSm·ωe,j·De,j(x,y,z) (12)
where x, y, z represent position coordinates, e represents the energy layer index, j represents the pencil beam spot index, SmIs a sensitive factor, omega, of the monitoring system related to the beam energye,jIs the weight of the jth beam of energy layer e, De,j(x, y, z) is the dose of energy layer e jth beam vs. position (x, y, z).
(2) Constructing a range uncertainty model of proton heavy ion intensity modulated radiotherapy
In this embodiment, the equations (10), (11), (13), (14), (15), (16), (17) described represent the uncertainty model equations and the equations used by the optimization algorithm, where the derivation of (12), (13) and (14) of equations (10) and (11) and the derivation of (15), (16) and (17) are independent;Represents a protonHeavy ion beam set, b (i) pencil beam set representing proton heavy ion beam ithen all pencil beams are countedRange of j pencil beam of ith proton heavy ion beamIt is shown that,indicating the normal range of the gas flow in the gas flow path,Indicating a maximum above/below normal range. The range of the jth pencil beam of the ith beam is expressed asFor arbitrary beamby usingrepresenting the set of possible ranges after uncertainty dispersion of pencil beam j e b (i). Pencil beams j andeach element in the array is in one-to-one correspondence, so that the array can be at any rangeby indexingand (6) unique confirmation. Use ofTo representRange index elementThe vector of (2). For convenience of definitionis at normal rangeIs used to determine the index of (1).
By usingVoxels representing the Target region (Target) and the Organ At Risk (OAR), respectively. By wjIndicating the intensity, theta, of the jth pencil beamvRepresenting the desired dose value for the ν th voxel.Is the j-th pencil beam (j e b (i)) at rangeDose contribution to the ν th voxel. The ν th voxel is subjected to a total dose of
Can be written as
after each pencil beam range introduces uncertainty, the objective function and constraints must depend on the set of possible ranges. The set of possible ranges can be considered as a uniform distribution, with the uncertainty model as follows:
Wherein R is the set of possible rangesthe number of range samples in the middle of the range,Andis a variable wjAt set of possible rangeThe k-th and p-th ranges,Andare respectively the rangeAndthe corresponding objective function value.
(3) Adaptive particle swarm robustness optimization method based on filtering algorithm
The method adopts a conventional particle swarm optimization process, namely particle initialization, fitness calculation, position updating and variation. And a new constraint processing method and a dynamic variation strategy are adopted for key links. Using filtering algorithms to "filter" solutions violating constraintsAnd (6) processing. In addition, in order to avoid the particle swarm optimization from falling into local optimization in the optimization problems of high variable dimension and complex form, an effective dynamic variation strategy is adopted. The velocity of the particle changes dynamically with the number of iterations, while the position of the particle changes with the change in velocity. The diversity function D and the fitness variance function σ are introduced below2。
f＝max{1,max{|f(i)-favg|}} 1≤i≤M
Wherein M is the population number and L is the search radius. f (i) is the fitness function value of the ith particle, favgAnd the average value of the fitness function values of the current particle swarm is obtained.
and respectively setting a threshold value for the population fitness variance function and the population diversity function, and disturbing the speed of individual particles of the population once the population fitness variance function value or the population diversity function value is lower than the threshold value continuously appears in the iterative process of the population. However, in the later stage of algorithm iteration, if the variation range is too large, the particles are greatly disturbed, and the search is influenced. Therefore, the particle velocity disturbance range needs to be continuously reduced, and the disturbance formula is as follows:
Wherein v isiis the velocity of the ith particle, t is the number of iterations, GmaxFor the total number of iterations, θ is a random number that follows a standard normal distribution N (0, 1).
(4) Systematic testing
the effectiveness of the algorithm is verified by 5 cases of the benchmark example, the proton intensity modulated radiotherapy clinical example of the head and neck tumor and the chest and abdomen tumor. The tumor radiotherapy group baseline geometric anatomical model (shown in fig. 2) was used, which is a standard model for testing the efficacy of different treatment regimens. The OAR is a circle with a radius of 1.5 cm, centered on a circular mold with a radius of 10 cm. The Tumor (Tumor) is shaped like a "horseshoe" and wraps half of the OAR, with an inner diameter of 1.8 cm and an outer diameter of 4 cm. The voxel resolution was 1 × 1 mm, with a total of 2186 Tumor voxels, 441 OAR voxels. The model is a uniform phantom of equivalent water. The proton beam was irradiated from three directions, 315 °, 0 °, and 45 ° from left to right, respectively. The beam spots were spaced 4 mm apart in both the axial (parallel) and transverse (perpendicular) directions for each proton beam, traversing the entire tumor and its extended 1 cm (isotropic) voxels. There are 877 beam spot positions in total for the three proton beams (different energy proton beams may have the same position).
the preferred embodiments of the invention disclosed above are intended to be illustrative only. The preferred embodiments are not intended to be exhaustive or to limit the invention to the precise embodiments disclosed. Obviously, many modifications and variations are possible in light of the above teaching. The embodiments were chosen and described in order to best explain the principles of the invention and the practical application, to thereby enable others skilled in the art to best utilize the invention. The invention is limited only by the claims and their full scope and equivalents.
Claims (3)
1. The proton heavy ion modulated intensity radiotherapy robust optimization method aiming at range uncertainty is characterized by comprising the following specific steps of:
Step 1: constructing a range uncertainty mathematical model based on the combination of optimality and robustness;
Step 2: adopting adaptive particle swarm robust optimization based on a filtering algorithm to obtain a robust optimal solution;
In step 1, the model embodies the optimality and robustness of radiotherapy plans caused by range uncertainty, and the specific formula is as follows:
In the above-mentioned formula,target function values under different ranges are represented, and optimality is embodied;The robustness is embodied, namely the deviation of target function values under different ranges; w represents the intensity of the pencil beam of the proton heavy ion beam, z represents the organ sampling point at rangeThe actual dosage to be subjected to the treatment,Representing a set of proton heavy ion beams, B (i) representing a set of pencil beams of beam i (i ∈ B), of the jth pencil beam of the ith beamfor firing rangeIt is shown that,indicating the normal range of the gas flow in the gas flow path,Indicating a maximum above/below normal range, the range of the jth pencil beam of the ith beam is indicated as For an arbitrary beam i ∈ B, use is made ofrepresenting the possible range set after uncertainty dispersion of the pencil beam j ∈ b (i), wherein R represents the number of elements in the range set; pencil beams j andEach element in the array is in one-to-one correspondence and has any rangeBy indexingunique confirmation; use ofindex element representing rangethe vector of (a); by usingRespectively representing a Target region (Target) and an Organ At Risk (OAR), v representing the v th voxel;AndWeights representing the target and the organs at risk, respectively; by wjIndicating the intensity, theta, of the jth pencil beamvA desired dose value representing a ν th voxel;
the step 2 comprises the following specific operations:
Step 2.1: firstly, carrying out filtering processing on a solution violating the constraint by adopting a filtering algorithm; the solution in formula (1) represents the weight of the pencil beam, the weight is a non-negative number, and if the weight obtained in the optimization process is negative, a filtering mode of setting zero is adopted for processing;
Step 2.2: an effective dynamic variation strategy is adopted, and a diversity function D and a fitness variance function sigma are introduced2；
f＝max{1，max{|f(i)-favg|}}1≤i≤M
wherein M is the population number, L is the search radius, f (i) is the fitness function value of the ith particleavgThe average value of the fitness function values of the current particle swarm is taken;
The variation function of particle velocity is shown below:
Wherein v isiIs the velocity of the ith particle, t is the number of iterations, Gmaxθ is a random number obeying a standard normal distribution N (0,1) for the total number of iterations;
Step 2.3: optimizing formula (1) by adopting a particle swarm algorithm, filtering solution violating the constraint by adopting the method in step 2.1 in the optimization process, and calculating the particles by adopting formula (2) and formula (3) to obtain D and sigma2If D and σ2if the threshold value exceeds the preset threshold value, performing mutation operation by adopting a formula (4); and after the stopping criterion is met, finally obtaining a robust optimal plan of the proton heavy ion intensity modulated radiotherapy aiming at the range uncertainty.
2. the robust optimization method of proton heavy ion modulated radiotherapy for range uncertainty of claim 1, further comprising the following sequence of operations before step 1:
1) Importing radiotherapy plan information
Radiotherapy planning information comprises CT images of a patient, beam parameters, accelerator hardware parameters, organ constraint parameters and range uncertainty parameters; a standard DICOM RT format is used as an input and output interface, a radiotherapy plan file of a commercial plan system is imported, after a radiotherapy plan is finished, proton heavy ion scanning point weight and dose information are exported to the DICOM file, and the supported standard files mainly comprise: CT, DICOM RT-Structure, DICOM RT-Plan, DICOM RT-Dose;
2) Calculating a dose contribution matrix
Calculating the dose contribution of proton heavy ion scanning points with different energies and different positions to all organ sampling points, adopting a fast analysis method-pencil beam dose calculation algorithm, and performing convolution superposition on the flux weight of the proton heavy ion at the position to be calculated and IDD data corresponding to the water equivalent depth, wherein the formula can be expressed as follows:
Dose(r，z)＝∑kΦk(r，z)·IDD(z) (10)
Wherein z represents the depth position along the beam, r represents the transverse position, k represents the index of the proton heavy ion pencil-shaped beam spot, Dose represents the total Dose in a voxel corresponding to a certain position, phi is the flux weight of the proton heavy ion, and IDD is the integral depth Dose data of the proton heavy ion pencil-shaped beam in the depth direction in water which is calculated in advance through actual measurement or a Monte Carlo method.
3. The robust optimization method of proton heavy ion modulated radiotherapy for range uncertainty as claimed in claim 2, wherein the specific steps of the sequence step 2) are:
Firstly, distributing points according to the full width at half maximum of a proton heavy ion beam spot in the beam direction; dividing each beam spot into 19 sub-beams; then tracking each sub-beam by adopting an improved Siddon algorithm, calculating the water equivalent depth corresponding to each voxel on a ray path, and looking up a table to obtain a corresponding IDD value; according to the Fermi-Eyges theory, a two-dimensional calibration mode is adopted to calculate the Gaussian distribution standard deviation sigma of the Torontal scattering corresponding to the voxel passed by the rayMCS(ii) a Then, the influence of nuclear scattering on beam flux transverse scattering is calculated according to a theoretical formula provided by Soukup and the like to obtain the weight omega of the nuclear scatteringNSAnd the standard deviation sigma of its Gaussian distributionNS(ii) a Considering the truncation error δ of the gaussian distribution, the total flux weight formula is obtained as:
Φ(x，y，z)＝(1-ωNS)·[G(x，y，σMCS)+δMCS]+ωNS·[G(x，y，σNS)+δNS] (11)
according to beam monitoring system sensitivity factor SmAnd superposing all the results according to weight to obtain the final dose (x, y, z):
dose(x，y，z)＝∑e∑jSm·ωe，j·De，j(x，y，z) (12)
where x, y, z represent position coordinates, e represents the energy layer index, j represents the pencil beam spot index, SmIs a sensitive factor, omega, of the monitoring system related to the beam energye，jis the jth energy layer eWeight of pencil beam, De，j(x, y, z) is the dose of energy layer e jth beam vs. position (x, y, z).
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710707207.6A CN107551411B (en) | 2017-08-17 | 2017-08-17 | Proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710707207.6A CN107551411B (en) | 2017-08-17 | 2017-08-17 | Proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107551411A CN107551411A (en) | 2018-01-09 |
| CN107551411B true CN107551411B (en) | 2019-12-17 |
Family
ID=60975711
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710707207.6A Active CN107551411B (en) | 2017-08-17 | 2017-08-17 | Proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107551411B (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3549636B1 (en) * | 2018-04-03 | 2021-03-24 | RaySearch Laboratories AB | System and method for radiotherapy treatment planning and delivery |
| EP3572124B1 (en) * | 2018-05-21 | 2024-05-08 | RaySearch Laboratories AB | System, computer program product and method for ion-based radiotherapy treatment planning |
| CN110570923B (en) * | 2018-06-06 | 2022-05-10 | 北京连心医疗科技有限公司 | Mixed Monte Carlo radiotherapy reverse optimization method, equipment and storage medium |
| CN111551978B (en) * | 2020-05-08 | 2022-09-27 | 中国辐射防护研究院 | Method for evaluating uncertainty of radioactive soil sample sampling process |
| CN114355768B (en) * | 2020-10-13 | 2023-12-15 | 天津工业大学 | Continuous annealing process robust optimization method based on improved multi-target dragonfly algorithm |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105787256B (en) * | 2016-02-03 | 2018-07-13 | 中国科学院合肥物质科学研究院 | A method of the particle beams fluence based on optimization algorithm is decomposed into pencil beam |
-
2017
- 2017-08-17 CN CN201710707207.6A patent/CN107551411B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN107551411A (en) | 2018-01-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107551411B (en) | Proton heavy ion intensity modulated radiotherapy robust optimization method aiming at range uncertainty | |
| US9844684B2 (en) | Optimization methods for radiation therapy planning | |
| Yepes et al. | A GPU implementation of a track-repeating algorithm for proton radiotherapy dose calculations | |
| CN102921115B (en) | Method for establishing measurement data-based simple and convenient irradiation source model of medical linear accelerator | |
| WO2008134869A1 (en) | Method and system for optimizing dose delivery of radiation | |
| CN110415785A (en) | The method and system of artificial intelligence guidance radiotherapy planning | |
| CN103282967A (en) | Automated treatment planning for radiation therapy | |
| CN108415058A (en) | The dose calculation methodology and system of radioactive ray | |
| CN107073286A (en) | Radiation therapy planning device, radiation therapy planning method and radiation treatment systems | |
| Sir et al. | Stochastic programming for off-line adaptive radiotherapy | |
| CN110310720A (en) | A kind of Monte Carto dosage computing method, equipment and storage medium | |
| Barragán Montero et al. | Performance of a hybrid Monte Carlo‐Pencil Beam dose algorithm for proton therapy inverse planning | |
| Schumann et al. | From prompt gamma distribution to dose: a novel approach combining an evolutionary algorithm and filtering based on Gaussian-powerlaw convolutions | |
| US20140275703A1 (en) | Method and system for dose determination of radiation therapy | |
| Zhang et al. | A two-stage sequential linear programming approach to IMRT dose optimization | |
| Moore et al. | Comparisons of treatment optimization directly incorporating random patient setup uncertainty with a margin‐based approach | |
| Winterhalter | Protons do play dice: validating, implementing and applying Monte Carlo techniques for proton therapy | |
| Zhong et al. | Generation of a novel phase‐space‐based cylindrical dose kernel for IMRT optimization | |
| Li et al. | Adaptive particle swarm optimizer for beam angle selection in radiotherapy planning | |
| Bogner et al. | Fast direct Monte Carlo optimization using the inverse kernel approach | |
| CN110289075A (en) | A kind of direct Ziye optimization method and system based on fuzzy entropy | |
| Bedford | A discrete ordinates Boltzmann solver for application to inverse planning of photons and protons | |
| Tian et al. | Moving GPU‐Open CL‐based Monte Carlo dose calculation toward clinical use: Automatic beam commissioning and source sampling for treatment plan dose calculation | |
| Chetty et al. | A Phase-Space Model for Simulating Arbitrary Intensity Distribution for shaped Radiosurgery Beams Using the Monte Carlo Method | |
| Siggel | Concepts for the efficient Monte Carlo-based treatment plan optimization in radiotherapy |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |