Phys. Med. Biol 43 (1998) 2109-2121 © 1998 IOP Publishing Ltd

Developing a dose-volume histogram
computation program for brachytherapy

E Panitsa1, J C Rosenwald2 and C Kappas1
1 Medical Physics Department, University of Patras, 26500 Patras, Greece
2 Service de Physique Médicale, Institut Curie, 26 rue d'Ulm, Paris 75005, France

Abstract.

A dose-volume histogram (DVH) computation program was developed for brachytherapy treatment planning in an attempt to benefit from the DVH's ability to present graphically information on 3D dose distributions. The program is incorporated into a planning system that utilizes a pair of orthogonal radiographs to localize the radiation sources. DVHS are calculated for the volume of tissue enclosed by an isodose surface (e.g. half the value of the reference isodose). The calculation algorithm is based on a non-uniform random sampling that gives a denser point distribution at the centre of the implants. Our program was tested and proved to be fast enough for clinical use and sufficiency accurate (i.e. computation time of 20 s and less than 2% relative error for one point source, for 100000 calculation points). Ibe accuracy improves when a larger calculation point number is used, but the computation time also increases proportionally. The DVH is presented in the form of a simple graph or table, or as Anderson's 'natural' DVH graph. The cumulative DVH tables can be used to extract a series of indexes characterizing the homogeneity and the dose levels of the distribution in the treatment volume and the surrounding tissues. If a reference plan is available, the DVH results can be assessed relative to the reference plan's DVH.

1. Introduction

Dose-volume histograms (DVHS) are currently used in many radiotherapy departments as an external radiotherapy treatment planning tool that facilitates the evaluation of three dimensional (3D) treatment plans. Their main contribution to treatment planning is their ability to summarize data of the 3D dose distribution into a dose-volume graph. In particular, DVHs represent the relation between dose and the volume receiving this dose in a volume of interest (i.e. tumour, sensitive organs). They permit us to overview 3D dose distributions and, consequently, to evaluate a 3D plan according to criteria of dose level, dose uniformity or possible undesirable hot and cold spots (Drzymala et all 1991).

Brachytherapy treatment planning can take advantage of evaluation tools such as DVHs, initially developed for external radiotherapy. The 3D volumetric dose distribution that DVHs make available can be of great importance to the improvement of brachytherapy planning (Wu et al 1988, Saw and Suntharalingam 1991). DVHs can be calculated for a specific anatomical structure such as the tumour or sensitive organs. Such an application of DVHs requires 3D structure delineation. This is only possible for CT-aided brachytherapy treatment planning systems ('1'PSs). DVHs have been calculated for such systems (Roy et al 1993, Ling et al 1994). They offer critical information on the dose levels in the target volume and on the occurrence of hot spots in sensitive organs.

However, many TPS's do not incorporate CT information. Instead, they use the classical planning method that utilizes a pair of radiographs to define the coordinates of the radiation sources. This method does not permit structure delineation. In those cases, the DVHs cannot be calculated for anatomical structures. DVHs are then calculated for the whole implant. It is then, more convenient to calculate the DVHs for the region enclosed by a specific isodose surface such as the reference isodose (i.e. the 60 Gy isodose according to ICRU 38 (1985)).

Developing a DVH computation program for brachytherapy is a procedure more or less based on the experience gained from the development of DVH computation algorithms for external radiotherapy treatment planning. These algorithms (Drzymala et al 1991) calculate the dose for a sufficiently large number of points in the volume of interest (i.e. anatomical structure or isodose contour). The calculation points are chosen through a sampling technique. Two types of sampling technique have appeared in the literature:

  • (i) Grid sampling, where the coordinates of the calculation points are increasing periodically, 50 that the points simulate a regular Cartesian grid.

  • (ii) Random sampling, where the point coordinates are chosen randomly.

An extensive discussion on the performance of the two techniques has already taken place (Dizymala et al 1991, Niemierko and Goitein 1990, 1993, Lu and Chin 1993a, b, Jackson et al 1993). The major conclusion stresses the fact that even though grid sampling is in general more accurate than random sampling when a sufficiently large number of calculation points is used, the last is preferred for its speed, reliability and flexibility. Random sampling can give acceptable results for a relatively small number of points (i.e. 400 quasi-random points are sufficient for an error smaller than 5% in external radiotherapy DVHs (Niemierko and Goitein 1990)), permitting high speed programs to be developed. It also gives uncorrelated point placement for ah geometries, offers the possibility to increase dynamically the number of calculation points and to have spatially non-uniform point distributions 50 as to capture high dose gradients (Kooy et al 1993).

The DVH representations in brachytherapy have taken two basic forms; either a graphic representation of a simple dose-volume graph, or of a dose-volume graph in an altered scale. Logarithmic scales have been used for brachytherapy DVHs, but the most convenient scale has been proposed by Anderson (1986) under the name of 'natural' DVHs. The 'natural' DVH has already been used for brachytherapy treatment planning evaluation (van't Riet et ai 1993, Low and Williamson 1995, Langmack and Thomas 1995).

Considering the above, we tried to develop and evaluate a DVH computation program for brachytherapy based on random sampling. Our program will be incorporated in a treatment planning system that utilizes the method of two orthogonal radiographs to localize the radiation sources. As there are no delineated structures, the DVH will be calculated for the region enclosed by a specific isodose surface. The DVH will be presented under the forms of graph, table, and 'natural' histogram.

2. Method and materials

The DVH is actually calculated as a two-column table representing dose values and the volume receiving these doses. In particular, two types of DVH are calculated: (i) the differential DVH, where the volume Vd (i) receiving dose in a dose interval, [D(i) - dD/2, D(i) + dD/2], is given as a function of the central dose value, D(i), of that interval, and (ii) the cumulative DVH, where the volume Vd(i) receiving dose higher than the lower limit of dose interval, [D(i) - dD/2, D(i) + dD/2], is given as the central dose value, D(i), of that interval.

The first column of the DVH table contains the central values, D(i), of the dose intervals for which the corresponding volumes will be calculated. The second column contains the volumes, V(i), that are found to receive doses in each dose interval. These volumes are actually found by calculating the dose on a number of points, Np, and converting the number of points, N(i), that receive doses in each dose interval to the volume, V(i), they represent.

The points where dose is calculated are chosen randomly in the volume of interest. A pseudo-random number generator is used to define a series of pseudo-random points which become the coordinates of the calculation points. The generator is a RANDU (IBM) generator that follows the formula Xn = 65 539Xn-1 MOD(232) to calculate the nth random it Xn from the (n-1)th random point Xn-1 where Xn-1 MOD(232) symbolizes the remainder of the division of Xn-1 by 232. The points that are generated with the above formula can be considered random since the period of the generator is much larger than number of points required for the DVH computation. In addition, different series of random numbers are generated at each run of the program, as the initializing number, X0, function of the computer clock.

The generated numbers are used to define the coordinates of the calculation points in a spherical coordinate system where the coordinates (r, q , j ) can take values from [0, 0, 0] to [R, 2p, p]. Spheric coordinates are used in order to produce a spatial non-uniforrn point distribution characterized by a denser distribution at the centre of the implant where the target volume usually lies. At the implant centre, high doses are accumulated and the dose gradient of each source is higher. The non-uniform point distribution has as a result the dependence of the volume that corresponds to each calculation point from the position of that point. Specifically, the volume DV that corresponds to a point (r, q , j ) is given by formula: DV = 2p2r2 sin j R/Np, where Np is the number of calculation points.

The DVH computation method should involve the dose calculation at a series of points lying in the space enclosed by two specific isodose surfaces. These isodose surfaces are defined by the user as the lower, Dmin and the upper, Dmax, limit of the dose axis of the DVH graph. However, the shape of the isodose surfaces is not known. Therefore, a sphere large enough to enclose both isodose surfaces is chosen as the volume in which the random points will be generated. Choosing a sphere as the sampling space complies with the use of spheric coordinate system. The sphere is actually defined by its centre, C, and its radius, The centre, C, is the geometrical centre of all sources, weighted by their individual a activities (case of linear sources) or by the square root of their individual activities (case of point sources). The origin of the coordinate system is placed at C. The radius, R, is the distance from C where the dose is definitely smaller than Dmin This condition is examined calculating the dose at points lying either at the extension of the line CPn, where Pn is most distant source from the centre (case of point sources), or at the extension of the geometrical axis of the most distant source (case of linear sources). Dose is expected to crease as we leave C. The point where the dose is found smaller than Dmin is the point that determines the radius, R, of the calculation sphere. A safety margin of 0.5 cm is used as to ensure that R is sufficiently large to include all points receiving dose larger than Dmax, no matter what the source configuration.

The data that the program requires as input are the following:

  • (i) The position of the radiation sources given either through a pair of radiographs or their coordinates in a Cartesian system.

  • (ii) The dose range (i.e. the lower, Dmin and upper, Dmax dose limit) for which the DVH will be displayed. Dmin is used to define the radius, R, of the sampling sphere as described above.

  • (iii) The number Ni of dose intervals (Ni < 1000). The width, dD, of the dose intervals can be easily extracted from the formula: dD = (Dmax - Dmin)/Ni.

  • (iv) The number of calculation points, Np. 50000 points are recommended as a lower limit.

The output of the program is the DVH (cumulative and/or differential) displayed as a graph or table, or as a 'natural' DVH.

Our program was developed in Fortran for VAX/VMS. The environment of the ISIS-1 treatment planning system was used. In particular, we used the source localization and dose calculation algorithms of ISIS-1 (Rosenwald 1975). Each source consists either of a point or a straight segment or a curve appr'0ximated by a series of straight segments. Its 'activity' is expressed in terms of reference air kerma rate. The dose at any point is calculated on the basis of the inverse square law with appropriate correction for attenuation by the tissue, assumed to be water equivalent (Pierquin and Marinello 1992). There is provision for automatic identification and reconstruction of curved wires.

3. Results

Several tests were made in order to check the accuracy and flexibility of the DVH computation program. Cases of one point source were firstly examined. These cases give the possibility to account for the accuracy of the program since it is relatively easy to extract manually theoretical results and compare the DVH to them. The theoretical calculations of the dose rate D' were executed manually using the equation
                 D'(r) = S/r2             (1)
where
hdvbf1.gif (10564 octets)       (2)

with k'R the reference air kerma rate (in µGy m2 h-1) (CFMRI 1983, Nath et al 1995), hdvbf2.gif (5071 octets)the mass energy absorption coefficient of water relative to that of air, and j (r) the effective transmission factor of water at distance r (in cm) from the source. j (r) is given by the formula of Meisberger et al (1968). The above formulae are followed by the treatment planning system in which our DVH computation program was incorporated.

Point sources of 192Ir of various reference air kerma rates were chosen in order to test the program. Taking into account that in the case of one point source the isodose curves take the shape of a sphere, the distance, r, can be related to the volume, V, of tissue receiving dose rates equal or greater than D'(r) by the following equation:

V = 4/3pr3 (3)

Combining equations (1) and (3) gives the cumulative DVH of one point source:

V = 4/3pS3/2D'(r)-3/2 (4)

Taking the theoretical values as reference, the errors of the program results were calculated. Several point sources were examined, in order to verify the accuracy of the program for both low and high dose, D, and distance, r, values. In figures 1(a) and 1(b), the calculated errors are presented for one point 192Ir source of reference air kerma rate 500 µGy h-1 m2, treatment duration of 1 hour and for a variety of calculation point numbers. Figure 1(a) presents the errors for low doses (i.e. far from the source, 2.6 cm < r <9.3 cm) while figure 1(b) presents the errors for high doses (i.e. close to the source, 0.4 cm < r <3.4 cm).

hdvbi1.gif (197637 octets)

The performance of our DVH computation program was also examined for some typical 'Paris system' source configurations. As an example, figure 2(a) presents a two plane 6 x 3 cm2 implant consisting of ten linear 192Ir sources placed parallel to each other. The sources have an active length of 3 cm and a linear reference air kerma rate of 50 µGy h-1 m2 cm-1. The treatment duration was selected to be 14 hours which gave 60 Gy at 85% of the basal dose of the configuration. The cumulative DVH for the configuration of figure 2(a) is presented in figure 2(b) as a graph and as a table. In these tests, it was not possible to use theoretical calculations as reference. However, repeating the calculations for a variety of calculation points showed that the results converge as the calculation point number increases. In particular, the results stabilize (variations less than 1%) at the 500000 point case. Therefore, the 500000 point case was selected as a reference for studying the characteristics of the computation for smaller calculation point numbers. The dose-volume table for a varying calculation point number was calculated for the configuration of figure 2(a). As an indication of the accuracy of the results, the relative variations comparing to the volume results for 500000 calculation points was calculated. The variations are presented in figure 2(c).

hdvbi2.gif (180990 octets)

4. Discussion

Summarizing the tests, we have to mention that the relative error of the program results as compared to theoretical calculations for one point sources is sufficiently low (i.e. less than 3%) for 50000 calculation points. The errors decrease as the number of the calculation points increases, being less than 1% for 500000 calculation points (figures 1(a), (b)). The program performs the same for both low and high doses, which ensures the reliability of the computation at all possible distances from the radiation sources. The one point source cases are characterized by high dose gradients; therefore, their performance is crucial for the assessment of the program accuracy.

Additional tests were performed for more complex source configurations. It was not possible to use theoretical results to test these configurations, but the 500000 points proved to be critical for the stabilization of the results (i.e. variations less than 1%). Using the 500000 point case as reference, it was found that the 50000 point cases maintain the variation relatively to the reference case below 3% (figure 2(c)). The number 50000 was chosen as the lower limit of the calculation point number, Np, accepted by the program.

The DVH computation program actually calculates DVH for the volume enclosed between two isodose surfaces defined by the user. The external isodose surface is defined by the lower dose limit, Dmin. This is usually a relatively low dose which permits the evaluation of both the treatment volume and the surrounding healthy tissue.

The DVH is displayed in the form of a graph or a table. The graph representation is not appropriate for clinical use because the predominance of the inverse square law overwhelms any other DVH characteristic. DVHs of configurations presenting the most different dose distributions are surprisingly similar. Therefore, the display of the actual dose-volume table is chosen as a more reliable representation form.

Important characteristics of the distribution can be extracted from the cumulative DVH tables. Such characteristics are the following: (i) the treatment volume, Vtr, which is the volume enclosed by the reference isodose, (ii) the variation, DVa, of the treatment volume from the treatment volume of a reference implant 'a', (iii) the fraction of volume that receives doses between 100% and 150% of the reference dose (i.e. dose homogeneity index, DHI) which is characteristic of the homogeneity of the distribution in the treatment volume, (iv) the fraction of volume receiving doses between 50% and 100% of the reference dose (i.e. health tissue dose index, I{m')I) which is characteristic of the dose distribution at the healthy tissue surrounding the implant and (v) the fraction of volume receiving doses higher than 200% of the reference dose (i.e. overdose index, ODI) which is characteristic of the excess of high dose in the treatment area. The above fractions are similar to concepts based on target volumes that have already been used for brachytherapy evaluation (Saw and Suntharalingam 1991). In our method, all indexes are related to the treatment volume (i.e. the volume of the reference isodose) and not to the target volume which cannot be determined. Therefore, the DVH program can only be useful in compliance with the isodose distribution and only if it has been verified that the reference isodose encloses the target.

A series of implants were assessed using DVH tables. The implants are rectangles 7 x 10 cm2 that follow the 'Paris system' rules. Their individual characteristics are presented in table 1. All implants have equal source strengths except for implants 'e'-'g' which have one of their sources with higher or lower strength. Implant 'c' has three planes of sources while all others are two plane implants. Implant 'd' differs from implant 'a' by small displacements of the sources simulating clinical implants. The isodose distributions (i e. 75, 50, 40, 30, 20 and 10 Gy) at the x-y plane are presented for all implants in figure 3. The target volume is not delineated. However, all implants aim to treat a volume of thickness at least 4 cm along the central line x = 0, width 7 cm and length 10 cm. The dwell times were chosen 50 that the 50 Gy isodose enclosed this volume. The only exception is implant 'g' where the dwell time is higher than necessary in order for the 50 Gy isodose to cover the portion of tissue close to the low strength source which is not covered in implant 'f' (i.e. implant 'f' is similar to 'g' but its dwell time results in an isodose contour of 50 Gy that counts for at least 4 cm thickness, 7 cm width and 10 cm length but does not enclose a significant portion of tissue close to the low strength source). DVH tables were calculated for the implants of table 1. The DVH covered a dose range from 5 to 130 Gy with 25 dose intervals and 50000 calculation points. Table 2 presents the indexes Vtr, DVa, DHI, HTDI and ODI that were extracted from the cumulative DVH tables for the implants of table 1. Bearing in mind that at all implants (except 'f') the 50 Gy isodose encloses the target, we can compare the dose distributions by means of dose levels and dose homogeneity using the results of table 2. As an example, we can conclude that since implant 'c' gives higher DHI, lower ODI and Vtr, and higher HTDI than implant 'a', it presents better homogeneity in the treatment area, gives a smaller overdose expose, irradiates a smaller volume to the reference dose and it also gives a higher fraction of dose to the surrounding healthy tissues.

hdvbi3.gif (196230 octets)

Another way of evaluating a brachytherapy plan is to compare its DVH table to DVH tables of reference cases (e.g. a theoretical implant is used as the reference case for examining the characteristics of an irregular implant applied clinically). Such comparison could show the dose levels in the implant, the presence of bot and/or cold spots and the relative homogeneity of the distribution in the treatment area. Figure 4 presents the variations calculated by the DVH tables of the implants of table 1 relative to the volume results of the reference implant (i.e. implant 'a' of table 1). These figures can be used to assess the dose distribution of the implants of table 1. As an example, implant 'c' seems to treat a smaller volume with the 50 Gy isodose than implant 'a'. It is also characterized by lower levels of dose in the treatment volume which results in better homogeneity and lower overdose expose, while it gives rather higher levels of doses to the healthy tissues surrounding the treatment volume.

hdvbi3c.gif (161509 octets)

An alternative representation form that overcomes the problem of the predominance of the inverse square law at the brachytherapy DVH has been proposed by Anderson (1986) under the name of 'natural' DVH. In the 'natural' DVH, an altered scale (i.e. -3/2 power of dose) is chosen 50 as to eliminate the inverse square law effects for the case of one point source. The 'natural' DVH is integrated in our program as a form of presentation. Several tests were done in order to check the validity of the 'natural' histogram representation when random sampling is used. Our results verified the accuracy of the calculations; however, the 'natural' graphs were in some cases unsuitable for clinical interpretation, due to great statistical variations. The need of holding the statistical variation below 5% became evident. More careful examination of the behavior of the variation showed that it is independent of the dose and of the number of sources but it depends on the number of calculation points per dose interval. We performed some tests using 40 dose intervals and a gradually increasing calculation point number in order to estimate the number of points per dose interval that is necessary to keep the statistical variation smaller than 5%. The test were performed with a source of 75 µGy h-1m2, a treatment duration of 1 hour and a dose range of 0.3 to 3 Gy. The dependence of the statistical variation on the number of points per interval is presented in figure 5. A number of 10000 points per interval was found necessary to keep the variation below 5%. This number is of the same order as in Anderson examples (1986), where grid sampling was used.

hdvbi4.gif (26675 octets)

Niemierko and Goitein (1994) proposed a DVH computation method where the histograms, named dose-volume distributions (DVDs), are calculated by sorting all the points according to their dose values. This method was also applied to our program as an alternative computation method. The calculated cumulative and differential DVD for one point source cases were found to be in complete agreement to the results of the traditional DVH.

hdvbt1.gif (156665 octets)

As expected, the running time of the DVH computation program was found to be proportional to the number of calculation points. As an example, the running time of the computation part of the program was found to be 20 s for 100000 points, for a Vax6000. This time slightly increases (by 10 s) when 'natural' DVHs were calculated. In contrast, the execution time of the DVD method depends mainly on the sorting algorithm. For an N2 sorting algorithm the running time is relatively long (i.e. 200 s for 1000 points) and proportional to the square of the calculation points. A serious improvement appears when an N lnN (i.e. 'quicksort') sorting algorithm is used. The execution time is then of the same order as in the traditional DVH computation method.

hdvbi5.gif (24512 octets)

5. Conclusion

A DVH computation program devoted to brachytherapy was developed in order to be used as a plan evaluation tool. The program was incorporated into a brachytherapy treatment planning system which utilizes two orthogonal radiographs to localize the radiation sources. Random point sampling was selected to keep the computation time low without losing in reliability. A spherical random point distribution was chosen in order to gain accuracy close to the centre of the implants where the treatment volume lies. The DVH is calculated for the volume of tissue enclosed by an isodose surface (i.e. if the reference isodose is applied then the DVH of the treatment volume is acquired). The program is accurate and fast enough for clinical use. In particular, the accuracy of the program is improved as the calculation point number increases but the computation time also increases proportionally. As an example, 100000 calculation points require 20 s of computation time and give error less than 2% for a one point source configuration.

The DVHs (cumulative and/or differential) are presented either in the form of a graph or as a dose-volume table. The table presentation is used in an attempt to overcome the predominance of the inverse square law on the results. It permits us to extract characteristics of the dose distribution that account for the homogeneity and the dose levels in the treatment volume and the surrounding tissues. It also permits the comparison to reference plans. Natural DVHs are also displayed but a large number of calculation points per interval is required in order to keep the statistical variation low (i.e. 10000 points per interval result in variation of less than 5%).

References
  1. Anderson L L 1986 A 'natural' dose volume histogram for brachytherapy Med Phys. 13 89~903
  2. Comite Francais pour la Mesure des Rayonnements Ionisants (CFMRI) 1983 Recommendations pour la deternnination des doses absorbées en curietherapie Rapport CFMRI 1 (Paris: Bureau National de Metrologie)
  3. Drzymala R E. Mohan R, Boewster L. Chu J, Goitein M, Harms W and Urie M 1991 Dose-volume histograms Int. J Radiat. Oncol Biol Phys. 21 714
  4. International Commission on Radiation Units and Measurements (ICRU) 1985 Dose and volume specification for operating intracavitary therapy in genecology ICRU Report 38
  5. Jackson A. Mohan R and Baldwin B 1993 Letters to the editor on 'Sampling techniques for the evaluation of treatment plans' Med Phys. 20 I 375-6
  6. Kooy H M, Nedri L A. Alexander E 3rd, Loelller J S and leedoux R J 1993 Dose volume histogram computations for small intracranial volumes Med Phys. 20 755-60
  7. Langmack K A and Thomas S J 1995 The application of dose-volume histograms to the Paris and Manchester systems of brachytherapy dosimetry Br. J Radiat. 65 42-8
  8. Ling C C, Roy J N, Sahoo N, Wallner K E and Anderson L L 1994 Quantifying the effect of dose inhomogeneity in brachytherapy: application to permanent prostatic implant with 125I seeds InI. J. Radiat Oncol Biol Phys.28 971-8
  9. Low D A and Williannson J F 1995 The evaluation of optimized implants for idealized implant geometries Med Phys. 22 1477-45
  10. Lu X-Q and Chin L M 1993a Sampling techniques for the evaluation of treatment plans Med Phys. 20 I5I~1 -1993b Letters to the editor on 'Sampling techniques for the evaluation of treatment plans' Med Phys. 20 1381-5
  11. Meisherger L L. Keller R J and Shalek R J 1968 The effective attenuation in water of the gamma rays of gold 198, Iridium 192, caesium 137, radium 226, and cobalt 60 Radiology 90 953-7
  12. Nath R, Anderson L L. Luxton G, Weaver K A, Williamson J F and Meigoeni A S 1995 Dosimetry of interstitial brachytherapy sources: recommendations of the AAPM Radiation Therapy Committee Task Group No 43 Med Phys. 22 209-34
  13. Niemierko A and Goitein M 1990 Random sampling for evaluating treatment plans Med PhyL 17 753-62 -1993 Letters to the editor on 'Sampling techniques for the evaluation of treatment plans' Med Phys. 20 1377-80 -1994 Dose-volume distributions: a new approach to dose-volume histograms in three-dimensional treatment planning Med Phys. 21 3-11
  14. Pierquin B and Marinello G 1992 Manuel Pratique de Curietlcrapie (Paris: Herman) pp 29-52 Rosenwald J C 1975 Automatic localization of curved wires used in brachytherapy-the COUREP program Comput. Programs Biomed 103-I2
  15. Roy J N, Walloer K E, Harrington P J, Ling C C and Anderson L L 1993 A CT-based evaluation method for permanent implants: application to prostate Int. J. Radiat. Oncol Biol Phys. 26 l63-9
  16. Saw C B and Suntharalingam N 1991 Quantitative assessment of interstitial implants Inl. J. Radiat Oncol Biol Phys. 20 135-9
  17. van't Riet A, te Loo H J, Mak A C, Veen R E, Ypma A F, Bos J, Hoekstra C J and Stenfert Kroese M C 1993 Evaluation of brachytherapy implants using the 'natural' volume-dose histogram Radiother. Oncol 26 82-4
  18. Wu A, Ulin K and Sternick E S 1988 A dose homogeneity index for evaluating 192Ir interstitial breast implants Med Phys. 15 104-7