Dose computation using Voronoi and Delaunay graphs

Michael Lahanas

Panitsa et al. used a non uniform distribution of sampling points for the calculation of DVH in brachytherapy. Sampling points are generated according to the following method:

Generate the spherical coordinated ( r, q, j) as three uniform random numbers within the intervals r Î [0, R], q Î [0, p] and j Î [0,2p]. The cartesian coordinates (x,y,z) of this point are given by x = rsin(q)cos(j), y = r sin(q)cos(j) and z = rcos(q).

The coordinate center is at the center of the source positions. The sampling point density decreases radial from a center simulating the falloff of the dose of a single source. The volume DV which corresponds to each sampling point is not constant but is a function of DV(r, q, j). The Jacobian determinant for the spherical coordinates transformation gives r2sinq. The angular and radial interval inside a sphere of radius R is 1/2p2R. It follows that the corresponding volume element DV for a sampling point, if Np sampling points are produced inside the sphere is: DV = 2 R p2r2sinq /Np

Unfortunately this method does not reduce significant the number of sampling points for the accurate calculation of DVHs.

Our stratified sampling method which uses partial uniform distributed sampling points, if large implants have to be considered, reduces the number of sampling points by a factor of 3-10. We expect to have the best result for non-uniformly distributed points. Both methods which we studied are based on the Voronoi diagram which is one of the most fundamental and useful constructs defined by irregular lattices.

We have developed methods to calculate DVHs form dose points which are non uniform distributed. We discuss the problem to reconstruct DVHs from these non uniform distributed dose points. In the last decade a progress was made in the non-uniform sampling theory by Feichtinger et al. who describe constructive and iterative algorithms to recover a band-limited function of d variables u with known spectral support from a set of samples which are sufficiently dense in the Euclidean d-space (compared to the size of the spectrum). The methods are all based on iterations, which start by producing first some kind of approximation to the unknown band-limited signal f in terms of certain auxiliary functions or measures. This can be a step function obtained by nearest neighborhood interpolation ( Voronoi method ) or piece-wise linear interpolation. However it can also be a discrete measure, a sum of Dirac measures sitting at the sampling positions. If one just takes the sampling values of f as amplitudes (maybe with a fixed relaxation parameter to compensate for loss of energy) one obtains the method suggested by Wiley-Marvasti, which is equivalent to the so-called frame method (the collection of shifted SINC-functions is a frame for the Hilbert space of band-limited functions with given fixed spectrum).

If one chooses the "volume" of the Voronoi-domain (half the distance between the two neighbors) as extra factor for each sampling value, one obtains a much better compensation for irregularities in the sampling set, in particular clusters are well treated by these adaptive weights. The method obtained in this way is called the ADP-method and turns out to be as efficient (in terms of iterations) as the Voronoi method, but much faster, i.e. cheaper in computational terms.

Meanwhile it is known, due to results of K. Gröchenig, that for the case of one dimensional distributions the classical Nyquist rate is also the critical maximal gap width, in the following sense: Whenever the maximal gap between subsequent sampling positions is strictly less than the Nyquist distance, at least the Voronoi and the adaptive weights method will be convergent, and even more, from the oversampling rate a prediction can be made about the rate of convergence of these algorithms. Even if some of these algorithms use special structures such as Toeplitz matrices the iterative nature of these methods is of limited practical use.

One approach would be to reconstruct the DVHs using directly the Voronoi diagram of the generators. This results in a tessellation of the space in regions called voronoi cells, which in the 3 dimensional case are polyhedral regions. If the density of the generators is locally sufficient high and fulfils the Nyquist criterion and if inside each such region the dose is sufficiently constant tne one approach to reconstructing the DVHs from the generators directly, could be by assigning to each Voronoi region a weight to the DVH which is proportional to its volume. In contrast in uniform point sampling the weight assigned to each sampling point is constant.

This approach is unfortunately computational complex, since it requires the determination of the Voronoi diagram not only as a set of vertices and edges which can be efficiently obtained using e.g. QHULL but also the construction of each Voronoi polyhedron. Additional the volume of each such polyhedron must be calculated. Even if this is simple an accurate determination using anatomical information requires that additional the part of the volume of each polyhedron inside each anatomical structure has to be known. Nevertheless using some approximate methods it remains to be seen how this method works.

An efficient way for the calculation of the dose-volume histograms is the use of sampling points adapted to the form of the dose distribution around each source. Since this is mainly given by 1/r2 the non-uniform distributed points must also be distributed as 1/r2 around each source.

For the production of the sampling points we use the following theorem:

If random numbers ri are uniformly distributed in [0...1] then si = F-1(ri) is distributed as F(x). We therefore produce radial distributed points around a source with n( r ) = 1/r2. The points are spherically symmetrical distributed in the range [rmin, rmax]. Then from:

We obtain the a point at radius r from the source where

Points on a sphere from these radial distributed points are produced using the following algorithm:

Choose a random number s uniformly distributed in [-1, 1] and a random number j uniformly distributed in [0, 2p]. The x, y and z coordinates of the point are given by:

The idea behind this algorithm is that for a sphere of radius r, the area of a zone of width h is always 2prh, regardless of where the sphere is sliced. Therefore the z-coordinates of random points on a sphere are distributed uniformly. Points obtained by Eq. 4 are then uniformly distributed on the surface of the sphere.

Using sampling points distributed according to equation 4 we have calculated dose volume histograms for a region which includes the cancer plus an additional margin, the so called PTV for an implant with approximately 30 sources.

The dose volume histograms obtained with the conventional method with uniformly distributed points is shown together with the histogram with 4836 Voronoi generators. The space is divided in Voronoi polyhedra in which the dose is assumed to be constant and given by the corresponding Voronoi generator. Only parts of the Voronoi polyhedra inside the volume where the dose volume histogram should be calculated is considered. A Monte Carlo integration routine is used to approximate this volume. The Voronoi polyhedra are calculated from the convex hull and the dual Delaunay diagram of the sampling points (i.e. Voronoi generators).

For the much larger body a region which includes and surrounds the PTV the dose volume histogram with the Voronoi method with 4819 generators describes even high dose ranges with high accuracy.

A test of the method can be obtained using the well known dose-volume histogram of a single source in a sperical region with radius r

The dose volume histogram is given for a single point like source of strength S by

dV/dD = -2pS3/2D-5/2

which follows from the relations

r = S1/2D-1/2

dr = -1/2S1/2D-3/2dD

dV = 4pr2dr = -2piS3/2D-5/2dD

We believe that the fluctuations in the DVH from the Voronoi method is due to undersampling. Using an interpolation we try to reconstruct know dose-volume histograms using nearest neighbor interpolation from the non uniform distributed generators. Goal of this approach is to minimize their number and to increase the number of interpolated points. The number of neighbors for each interpolated dose value remains is usual smaller than the number of sources. We expect to speed-up by this approach the dose calculation for large implants with many sources.

Some 3D Voronoi Cells with sampling points to calculate the volume inside the PTV

Sampling points used to generate the Voronoi diagram. Close to the sourcs the density increases.

Some of the Voronoi cells getting smaller close to the radioactive sources

The extreme complex Voronoi diagram in 3D in the PTV. The graph complexity is higher close to the sources where the dose gradients are also large.

and in the Body.. less dense as the dose decreases