Method, computer program product, and computer system for dose calculation

Abstract

A method for calculating a dose for brachytherapy is disclosed, in which radiation is provided by one or more radiation sources from at least a first dwell position and a second dwell position within a patient. The method includes performing a first dose calculation for the first dwell position according to a first dose calculation algorithm, performing a second dose calculation for the second dwell position according to a second dose calculation algorithm, and using the sum of the first dose calculation and the second dose calculation as a total calculated dose. By combining different dose engines with different characteristics, calculations with sufficient accuracy and good time efficiency can be provided.

Description

Technical Field

[0001] The present invention relates to radiotherapy treatment, and more specifically, to brachytherapy.

Background Art

[0002] Unlike external beam radiotherapy, where radiation is delivered to a patient from outside in the form of a beam, brachytherapy involves inserting a radiation source into or adjacent to the target volume within the patient. The radiation source can be in the form of a seed, ribbon, or capsule. Brachytherapy has several advantages, such as being able to place the radiation source close to the target, and the radiation source remaining in a fixed position relative to the target even if the patient moves as a whole or there is internal movement within the patient. Also, a high dose rate is possible, allowing for treatment in only a few fractions.

[0003] In some forms of brachytherapy, the radiation source is placed, moved, and inserted into the patient's body typically using one or more catheters by utilizing natural cavities or interstitial channels. The radiation source moves through the channel and stops at predetermined positions along a trajectory, known as dwell positions. The time spent at each dwell position, i.e., the dwell time, is determined individually for each position and can typically be on the order of 1 to 2 seconds.

[0004] As with any type of radiotherapy, it is important to calculate the dose delivered to different parts of the patient. Since the radiation source emits in all directions from every dwell position within the patient's body, a large amount of calculation is required. Generally, known algorithms for such calculations are either based on simplified assumptions that speed up the calculation but are not necessarily highly accurate, or are based on more accurate but time-consuming dose engine calculations.

[0005] One way to calculate the dose involves an analytical algorithm that uses the cylindrical symmetry of the dose distributed by the radiation source. This algorithm is based on the assumption that the surrounding tissue is equivalent to water and is homogeneous. This simplification means that the algorithm is fast, but it also means that it is not necessarily highly accurate, especially in tissues that are particularly inhomogeneous or have material properties that are very different from water. Relevant material properties include density, effective electron density, the photoelectric effect, and Compton scattering.

[0006] Another way to calculate the dose involves a Monte Carlo dose engine. Such dose calculations can be performed extremely accurately, but are currently extremely time-consuming and may not be suitable for clinical applications, especially with regard to optimizing treatment plans.

Summary of the Invention

[0007] The object of the present disclosure is to provide a dose calculation method for brachytherapy that is highly accurate and fast compared to prior art methods.

[0008] The present disclosure relates to a dose calculation method for brachytherapy in which radiation is provided by one or more radiation sources from at least a first dwell position and a second dwell position within a patient, and the radiation provided at each dwell position results in a dose contribution from the dwell position. The method includes selecting a first dose calculation algorithm used to calculate the dose contribution from the first dwell position and a second dose calculation algorithm, different from the first dose calculation algorithm, used to calculate the dose contribution for the second dwell position, calculating the dose contribution from the first dwell position according to the first dose calculation algorithm and calculating the dose contribution from the second dwell position according to the second dose calculation algorithm, and calculating the total dose as the sum of the dose contributions from at least the first dwell position and the second dwell position.

[0009] Such a method enables an optimal combination of dose calculation methods with respect to time and accuracy, taking into account the patient's anatomical structure. This method identifies positions that require more complex algorithms to obtain highly accurate results. A fast but less accurate algorithm can be used when such an algorithm returns sufficient results. Thus, the method saves time by using a faster method in regions where less accuracy is required, i.e., regions where the surrounding tissue has a density close to that of water and shows little variation, while returning highly accurate values in complex regions. Thus, sufficient results can be achieved with reasonable calculation times.

[0010] The first dose calculation algorithm can be an analytical algorithm that assumes a uniform density of the tissue surrounding the first dwell position, such as the density of water. Such algorithms include ray tracing and TG43, which are fast but not necessarily very accurate.

[0011] The second dose calculation algorithm is preferably configured to take into account inhomogeneities. Such algorithms are slower but yield more accurate results. Examples of such algorithms include MC dose engines and Boltzmann equation solvers.

[0012] The dose calculation algorithm for each dwell position is preferably selected based on an indicator value that indicates either the difference in material properties such as the density between water and the tissue surrounding the dwell position, or the variation in material properties in the tissue surrounding the dwell position, or a combination of these two. This facilitates the selection of an appropriate dose calculation algorithm for each dwell position.

[0013] The indicator value can be an inhomogeneity indicator for the dwell position, and the inhomogeneity indicator indicates the variation in tissue properties in the volume surrounding the dwell position. A faster algorithm is selected for regions with low inhomogeneity, and a slower and more accurate algorithm is used for regions with high inhomogeneity.

[0014] The non-uniformity index for each residence position can be calculated based on performing ray tracing in two or more directions near the residence position and determining the tissue characteristics of the intersection voxels for each ray tracing operation.

[0015] Alternatively, the non-uniformity index for each residence position is calculated based on defining a sub-volume containing a plurality of voxels around each residence position and determining the tissue characteristics in at least some of the plurality of voxels. Since the variation in tissue characteristics indicates the degree of non-uniformity, a large variation results in a high non-uniformity index, while a small or no variation results in a low non-uniformity index.

[0016] Instead of or in addition to the non-uniformity index described above, the dose calculation method for each residence position can be selected based on comparison with a uniform reference substance such as water. This makes it possible to select different dose calculation methods for different uniform tissues. Typically, tissues with a density different from that of the reference substance should use a more accurate dose calculation method even if the tissue has a low non-uniformity index.

[0017] Three or more dose calculation algorithms can be used. For example, high-speed, medium-speed, and low-speed but highly accurate algorithms can be used in regions of low importance, moderately important regions, and very important regions, respectively. Typically, extremely non-uniform regions are the most important regions, and extremely uniform regions can be processed using a high-speed algorithm.

[0018] The dose calculation method according to the present invention can be used in any suitable situation, but is particularly useful for optimization. Accordingly, the present disclosure also relates to an optimization method for brachytherapy for optimizing a set of residence times for a corresponding set of residence positions. The optimization method for brachytherapy includes calculating the total dose using the method according to any of the above embodiments and using the calculated total dose in an optimization procedure.

[0019] The present disclosure also relates to a computer program product configured to perform dose calculations according to a first dose calculation algorithm and a second dose calculation algorithm, the computer program product comprising computer-readable code means which, when executed within a computing device, cause the computing device to execute the method according to any of the above embodiments. As is common in the art, the code means may be stored in a non-transitory memory. The present disclosure also relates to a computer comprising a processor and a program memory, the program memory storing such a computer program product so as to be executable by the processor. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] The invention will be described in more detail below, by way of example, with reference to the accompanying drawings.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

[0021] FIG. 1 shows the head of a patient 10 undergoing brachytherapy in a side view, including the treatment site around a target 12 that is exposed to radiation. Three channels 14 are provided, through which a radiation source is passed in a manner known in the art. As will be appreciated, typically there are several channels, and each channel has several dwell positions. For clarity, only three channels are shown in this figure. The treatment plan includes several dwell positions within each channel 14 where the radiation source temporarily stops, and the dwell time for each dwell position indicating the length of time the radiation source stops at that position. The total dose is the sum of the doses delivered from all the dwell positions. To protect organs in a dangerous state, a shield (not shown) can be applied to stop radiation in a specific direction. As will be appreciated, the geometric shape of the treatment site is nearly uniform in some regions and more non-uniform in other regions due to air cavities, bone, teeth, and possibly other substances. This affects the dose delivered in different directions from radiation sources at different positions. The methods for doing this are known in the art and include determining the dwell positions and the nature of the radiation source. The geometric shape of the radiation source also affects the dose in different directions. Thus, the dose calculation algorithm for each dwell position can be selected according to the geometric shape of the treatment site near the dwell position and / or the geometric shape of the radiation source.

[0022] In some embodiments, a dose calculation algorithm that is fast but not very accurate is the commonly used analytical TG43 dose algorithm. Examples of dose calculation algorithms that are accurate but slow include dose engine algorithms such as Monte Carlo-based algorithms.

[0023] FIG. 2 is a flowchart of a method according to an embodiment of the present invention. The method is executed after identifying a target and any organs in a dangerous state (S21), and after identifying several dwell positions within the treatment area (S22). This constitutes the input data to the method.

[0024] The total dose obtained from the dwell positions is necessary during the optimization of the treatment plan, but also when calculating the final dose of the resulting treatment plan. The final dose is used to evaluate the prescription and clinical goals and to approve the delivery of the planned treatment. The accuracy requirements are higher for the final dose than for the dose used during optimization. According to the present disclosure, two or more different dose calculation algorithms with different characteristics are used. In step S23, an optimal dose calculation algorithm is selected for each dwell position. The selection of the dose calculation algorithm for each dwell position can be based on any suitable set of criteria. In some embodiments, the selection is based on an indicator value indicating a difference in material properties, such as the density between water and the tissue surrounding the dwell position, or a variation in material properties in the tissue surrounding the dwell position, or a combination of these two. The method for determining the indicator value will be described in more detail below.

[0025] In step S24, the dose contribution from the radiation source at each dwell position is calculated using the algorithm selected for that dwell position in step S23. In step 25, a dwell time is selected for each dwell position. In step S26, the total dose obtained is calculated as the weighted sum of all the dose contributions weighted by their respective dwell times at the dwell positions calculated in step S24 and obtained in step 25. The method can include calculating the dose contribution from dwell positions belonging to several channels or cavities within the patient.

[0026] The method according to FIG. 2 can be used as part of the optimization procedure, the dwell time is used as an optimization variable in dose-based optimization, and the dwell time is changed based on the realization of a user-specified optimization function that is equal across the total dose obtained. This is shown by an optional step S27 that allows the method to be repeated from step S25 until the optimization is complete, for example, when a predetermined number of iterations is reached, or when the optimization function is realized, for example, when an appropriate total dose is achieved.

[0027] Figure 3 shows a first method for calculating the above-described indicator values. Reference numeral 31 indicates the channel through which the line source passes. A plurality of stagnation points 33 are defined within the channel. For illustrative purposes, one stagnation point 33a has been selected. At this selected stagnation point 33a, several directions (three directions are shown in Figure 3) are defined or sampled, and ray tracing is performed in these directions. Ray tracing here includes accumulating terms for each intersection voxel i in a ray tracing that depends on the relevant tissue property, shown here as ρ. Suitable tissue properties to consider include tissue density, effective electron density, cross-sectional area of the photoelectric effect or Compton scattering. The resulting terms are labeled H1, H2, and H3 in Figure 3. Since the dose from the line source at each stagnation position decreases with increasing distance, this term preferably decreases with increasing distance from the stagnation position and can be the main cause for the material variation closer to the line source to have a greater impact on the dose distribution than the material variation far away. The indicator value for the selected stagnation point is calculated as a sum based on the results of these three ray traces, and in some cases as a weighted sum over all ray traces j. Since the stagnation positions are usually positioned close to each other, the indicator value can also be calculated for groups of adjacent stagnation positions rather than individually for each stagnation position.

[0028] The contribution to the indicator value from each ray trace in the example of Figure 3 can be expressed as Equation 1 below.

[0029]

Equation

[0030] The comprehensive index value H in this case tot can be expressed as Equation 2 below. H tot = Σ j H j (2) That is, it is the sum of the total index value contributions calculated according to Equation 1 for all directions j selected for the ray trace.

[0031] Alternatively, or in addition to the comparison with water, the comparison between the tissue properties of different voxels within the relevant volume can be used as the inhomogeneity index. In this case, Equation 1 is replaced by Equation 1a below.

[0032]

Number

[0033] If the tissue surrounding the residence position is relatively homogeneous, it may still be desirable to select a more accurate dose calculation method. Therefore, the selection of the dose engine can also be based on the difference in tissue density between the tissue surrounding the residence position and the reference density value, and a difference exceeding the set threshold indicates that a more accurate dose calculation method should be used. The appropriate reference density value is often the density of water.

[0034] FIG. 4 shows a second method for determining the index value H. Similar to FIG. 3, a channel 41 including several residence positions 43 is shown. In that case, the sub-volume 45 is defined around the residence position. One or more associated tissue characteristics are determined for all the voxels within the sub-volume, and the index value is calculated by accumulating the tissue characteristic values for these voxels. This procedure can be performed more efficiently by considering several adjacent residence positions together, i.e., by defining one sub-volume covering a plurality of residence positions, and / or by caching a part of the inhomogeneity index for each sub-volume. In this case, the comprehensive index value H tot is obtained by the following Equation 3.

[0035]

Equation

[0036] That is, the index value H tot is given as the sum over all the voxels within the sub-volume 45 of the difference between the tissue density of the voxel and the density of water divided by the square of the distance r between the voxel i and the center of the residence point. Similar to Equation (1), instead of density, several other substance characteristics, relative differences, or differences of a series of substance characteristics can be considered. Also, instead of the density of water, another uniform density can be used. Similar to the method including ray tracing, the comparison between the tissue characteristics of different voxels within the associated volume can be used as an inhomogeneity index instead of or in addition to the comparison with water.

[0037] Therefore, the selection in step S23 can be divided into sub-steps as shown in FIG. 5. In the first step S51, using ray tracing and Equations 1 and 2, or using the sub-volume and Equation 3, the index value H tot for the residence position or group of residence positions is calculated. In step S52, the index value H tot is compared with a threshold value. The index value H totIf it is lower than the threshold value, this indicates that the tissue characteristics are similar to those of water and / or are relatively uniform, meaning that a high-speed calculation algorithm can be used at the corresponding residence position even if it is not very accurate. The index value H tot If it is higher than the threshold value, this indicates that the tissue characteristics are different from those of water and / or are non-uniform, meaning that a more accurate but time-consuming calculation algorithm can be used at the corresponding residence position. As can be understood, a very high-speed algorithm can be used below the lowest threshold value, a medium-speed algorithm can be used between the lowest threshold value and a higher threshold value, and above the higher threshold value, several different algorithms and their respective threshold values can exist such that a slower but very accurate algorithm can be used.

[0038] In the simplest case, two calculation algorithms are used. One is a low-speed algorithm that is highly accurate when accuracy is required, and the other is a faster but less accurate algorithm for less important voxels. For the highly accurate calculation algorithm, a dose engine such as a Monte Carlo dose engine or a Boltzmann equation solver can be used. In the case of a fast but less accurate calculation algorithm, as described above, the analytical TG43 dose engine is most commonly used. The algorithms using the highly accurate dose engine and the less accurate dose engine can be combined according to their characteristics and the desired dose calculation speed. Examples of dose calculation algorithms that can be used are as follows. · Monte Carlo (MC). The accuracy of this algorithm is very high, but it is typically slower than other algorithms. The time increases or decreases depending on the size of the volume for which the dose is calculated and the number of approximations performed within the MC engine. The time is inversely proportional to the desired statistical uncertainty. · Boltzmann equation solver. The speed and accuracy are comparable to those of Monte Carlo. · Collapsed cone (CC). Medium speed. It can take into account the tissue inhomogeneity, but includes approximations that reduce the accuracy of the algorithm under certain conditions such as high inhomogeneity, high energy, or large treatment volumes. · Pencil beam. High speed, but low accuracy under inhomogeneous conditions. · SVD. Faster than pencil beam and CC. In terms of accuracy, it is comparable to the pencil beam. · Ray tracing. Very fast, but there are limitations to the accuracy under inhomogeneous conditions. · TG43. Currently the most commonly used. This algorithm is very fast, but tends to have reduced accuracy under inhomogeneous conditions or conditions different from water.

[0039] For each of the listed algorithms, different settings that affect performance and accuracy can be used. Combinations of different algorithms as described above can also include algorithms that are executed by the same dose engine but with different settings.

Claims

1. A computer-implemented dose calculation method for brachytherapy in which radiation is provided by one or more radiation sources from at least a first residence position and a second residence position within a patient, and the radiation provided at each residence position results in a dose contribution from the residence position, the method comprising: selecting a first dose calculation algorithm used to calculate the dose contribution from the first residence position and a second dose calculation algorithm, different from the first dose calculation algorithm, used to calculate the dose contribution for the second residence position; calculating the dose contribution from the first residence position according to the first dose calculation algorithm and calculating the dose contribution from the second residence position according to the second dose calculation algorithm; and calculating a total dose as the sum of the dose contributions from at least the first residence position and the second residence position.

2. The dose calculation method according to claim 1, wherein the first dose calculation algorithm is an analytical method assuming a uniform density of the tissue surrounding the first residence position, such as the density of water.

3. The dose calculation method according to claim 1 or 2, wherein the second dose calculation algorithm is configured to take into account the non-uniformity in the tissue surrounding the second residence position.

4. The dose calculation method according to any one of claims 1 to 3, wherein the dose calculation algorithm for each residence position is selected based on the geometric shape surrounding the residence position and / or the geometric shape of the radiation source.

5. The dose calculation method according to any one of claims 1 to 4, wherein the dose calculation algorithm used for each residence position is selected based on an index value indicating either a difference in material properties such as density between water and the tissue surrounding the residence position, or a variation in material properties in the tissue surrounding the residence position, or a combination of the two.

6. The method according to claim 5, wherein the index value for each residence position is based on a non-uniformity index of the residence position, and the non-uniformity index indicates a variation in tissue properties in the volume surrounding the residence position.

7. The method according to claim 6, wherein the non-uniformity index for each residence position is calculated based on performing ray tracing in two or more directions near the residence position and determining the tissue properties of the intersecting voxels for each ray tracing operation.

8. The non-uniformity index for each residence position is calculated based on determining a sub-volume including a plurality of voxels around each residence position and determining the tissue characteristics in at least some of the plurality of voxels. The method according to claim 6.

9. The dose calculation algorithm for each residence position is selected based on a comparison between the material characteristics of the tissue surrounding the residence position and a uniform reference material. The method according to any one of claims 1 to 8.

10. Performing a third dose calculation for a third residence position according to a third dose calculation algorithm, and using the sum of the first dose calculation, the second dose calculation, and the third dose calculation as the total calculated dose. The method according to any one of claims 1 to 9.

11. An optimization method for computer-implemented brachytherapy for optimizing a residence time set for a corresponding residence position set. The optimization method includes calculating the total dose using the method described in any one of claims 1 to 10, and using the total calculated dose in an optimization procedure. An optimization method for computer-implemented brachytherapy.

12. A computer program product configured to perform dose calculations according to a first dose calculation algorithm and a second dose calculation algorithm. When the computer program product is executed within a computing device, it comprises computer-readable code means for causing the computing device to execute the method according to any one of claims 1 to 11. A computer program product.

13. A computer program product comprising a non-transitory storage medium storing the computer program product according to claim 12.

14. A computer comprising a processor and a program memory, wherein the program memory stores the computer program product according to claim 12 so as to be executable by the processor.

Images