|
|
Michael Lahanas and Thorsten Kemmerer
Theory Preparation of the data input 1D case 3D case References Implementation
Theory Assume that we want to compute the dose distribution
This is the conventional approach. Is there another calculation method and possible faster? The underlying idea of computing the dose distribution D(x) for Ns sources with a FFT method is as follows. Given two functions f(t) and g(t), and their Fourier transforms F[f] and F[g] respectively, convolution f*g of the functions f and g is defined as
The convolution theorem states that the Fourier transform of the convolution is just the complex product of the individual Fourier transforms i.e.
For the one-dimensional case Eq. (1) becomes:
The discrete summation can be written as an integral:
where g(x)
describes the geometric distribution of sources and their associated strength in space described by the d(x) Dirac function. This Equation describes a convolution of the kernel f(x) with the function g(x) which describes the distribution of sources in space.
This convolution can be very efficiently calculated using FFT algorithms. From the convolution theorem it follows that the source distribution D(x) can be obtained by
where F-1 denotes the inverse Fourier transform and N is the sampling size. The extension in three dimensions gives for the three-dimesional dose distribution in space D(x,y,z)
where Nx,Ny,Nz denotes the number of sampling points in the x,y and z-direction respectively and which define the achieved resolution.
This Equation describes a method for the calculation of the three-dimensional dose distribution D(x,y,z) from the Fourier transform of the source distribution function g(x,y,z) and the Fourier transform of the dosimetric kernel f(x,y,z) of a single source. The dose distribution can be obtained from the inverse Fourier transform of their normalized complex product. In an optimization procedure this function changes due to variations of the strengths, while the kernel remains constant. With a discrete Fourier transform (FFT) it is therefore necessary to compute the FFT of f only once. The FFT of g, the product of the normalized FFTs of f and g and their inverse FFT has to be computed in each iteration step.
The time taken for the multiplication of F[f(x)] and F[g(x)] can be neglected, so that the computational time for the dose distribution is effectively the time taken for the calculation of the two Fourier transformation steps. While for the conventional method the dose distribution calculation time is proportional to the number of sources, for the FFT based convolution method it is independent of the number of sources.
Preparation of the data input
1D case
For discrete functions f and g with period N defined by their samples fi, gi, i=0,1,…,N-1, the discrete convolution f*g of length N is given by
The convolution is symmetric i.e. f*g = g*f. Two important issues have to be considered using the discrete convolution: the input format of the data and the convolution size.
The discrete convolution assumes that both functions are periodic and of the same length. The first assumption leads to wrap around effects which can be avoided by extension of these functions with zero padding. Since the FFT requires only positive indices, the application of Eq. (9) requires a shift of the function. From the assumed periodicity this can be achieved by a wrap around ordering of one of the two functions. Wrap around order means that the first half of the array contains the discrete function values at positive coordinates, while the second half of the array contains the function values at negative coordinates.
Figure 1 Schematic diagram illustrating the FFT based method for calculating the dose distribution for the one-dimensional case of three sources with different strengths. (a)-(d), (g) The f and g array are shown as used in the discrete convolution, with wrap around order and zero padding for avoiding overlaps at both ends.
The method is shown schematically for a one-dimensional example with three sources as shown by g(x) in Fig. 1(a) with the corresponding strength of each source. N sampling points are used for the function g(x) and extended on both sides with N/2 zeros resulting in an array of size 2N. For f(x) which is shown in Fig. 1(b) we use 2N sampling points. This procedure is necessary for the situation where sources are located close to one end of the computational grid and the dose contribution at a grid point at the opposite end of the array has to be calculated. f must therefore be sampled up to a distance equal to the size of the dose calculation grid.
In Fig. 1(c) the kernel f is shown in wrap around order necessary for the FFT. In Fig. 1(d) and 1(e) the real and imaginary parts of the Fourier transform of g and f respectively are shown. In Fig. 1(f) the real and imaginary components of the product of the Fourier transforms of Fig. 1(d) and Fig. 1(e) is shown. For the kernel f we use a 1/(r2+e2) form with e =0.01mm (see also Results section). The curvature in Fig. 1(e) of the real part of the Fourier transform of Fig. 1(c) is actually much smaller but a scale has been selected to make it visible. The inverse Fourier transform of this product is shown in Fig. 1(g). It is the dose distribution produced by the three sources of Fig 1(a) shown using a logarithmic scale. For comparison we include the dose distribution obtained by Eq. (1). The dose distribution D(x) can be obtained from N central points of the resulting array in Fig. 1(g) of size 2N. The abscissa of Fig. 1 is described as index in the range 1 to 2N. The ordinate of Fig. 1(b), 1(c) and 1(g) is dose rate in arbitrary units. Fig. 1(d) and 1(g) are very similar due to the very small curvature in Fig. 1(e) and the vanishing imaginary components of the real and symmetric kernel f.
3D case
For the three-dimensional extension of this method we have to consider the wrap around of the kernel f(x,y,z). The Fourier transformation assumes that the origin of the coordinate system is at the center of the data matrix, whereas the origin of the computational data array is at the lower left point show in Fig. 2(a). The three-dimensional extended array with 2Nx´2Ny´2Nz samples is filled by the discrete values of f. We need to rearrange the data in this array in wrap around order. In Fig. 2(a) eight partitions of the original f array are numbered to visualize the transformation of the data which is required for a correct implementation of the FFT based method. The eight partitions must be rearranged in the manner shown in Fig. 2(b). The array g is expanded similarly to a 2Nx´2Ny´2Nz array with zero padding but without rearrangement as shown in Fig. 3. The inverse FFT of the product of these two extended arrays results in an array of size 2Nx´2Ny´2Nz. The dose distribution is obtained from the Nx´Ny´Nz grid points centered inside this array.
Figure 2 Transformation of the function f(x,y,z) required for a correct application of the convolution based method for the calculation of the dose distribution in brachytherapy. (a) The original kernel f(x,y,z) with 2Nx´2Ny´2Nz sampling points shown with eight partitions. (b) The rearrangement of these partitions required for a correct application of the convolution.
Figure 3 The function g(x,y,z) with Nx´Ny´Nz sampling points. The convolution requires the use of an extended array of size 2Nx´2Ny´2Nz. The space outside the inner data block is filled with zeros.
References
A. L. Boyer, E. C. Mok, Brachytherapy seed dose distribution calculation employing the fast Fourier transform, Med. Phys 13, 525-529 (1986)
A. K. Erdi , E. D. Yorke, M. H. Loew, Y. E. Erdi, M. Sarfaraz, B. W. Wessels, Use of the fast Hartley transform for three-dimensional dose calculation in radionuclide therapy, Med. Phys. 25, 2226-2233 (1998)
J. W. Cooley and J. W. Tukey, An algorithm for the machine computation of the complex Fourier series, Mathematics of Computation” 19, 297-301 (1965)
C. Temperton, “A generalized prime factor FFT algorithm for any n=2p3q5r, SIAM Journal on Scientific and Statistical Computing 13, 82-87 (1992)
I. W. Selesnick and C. S. Burrus, Automatic generation of prime length FFT programs, IEEE Transactions on Signal Processing 44, 14-24 (1996)
J. L. Lawall, Faster Fourier transforms via automatic program specification, INRIA Report 3437 (1998)
M. Frigo and S. G. Johnson, The fastest Fourier transform in the west, Technical report of the MIT Laboratory for Computer Science MIT-LCS-TR-728 Cambridge, MA, (1997)
C. E. Shannon, A mathematical theory of communication, The Bell System Technical Journal 27, 379-423, (1948)
W. H. Press, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes in C: The Art of Scientific Computing, New York, NY: Cambridge University Press, 2nd edition. 1995
Implementation:
The file FFTConv.zip contains two Microsoft Visual C5.0 projects produced by Thorsten Kemmerer:
1) In the directory test 1D FFT the project 1D_mit_NR.dsw
The resulting dose distribution is written in the file out/ConvDosis.txt.
2) In the directory test 1D KON the project 1D_Konventionel.dsw
The resulting dose distribution is written in the file out/FFTDosis.txt.
We assume that a subdirectory exists at the executable directory
Project 1 contains an example of the conventional calculation of the dose distribution with 3 sources with arbitrary strengths 8, 4 and 6 respectively.
Project 2 contains the calculation of the dose distribution as in project 1 using 1D FFT. A numerical recipes FFT routine is used. Much faster is the Fastest FFT in the West library FFTW.
For the kernel f(x,y,z) we introduced a cut-off using a parameter e to avoid singularities. The effect of this parameter can be neglected as long as e is smaller than the grid constant.
References
PDF Files are not more available!!! Thesis Thorsten Kemmerer (in German) Consider that the corresponding figure 4.10 of Figure 2 of this document is not correct drawn in the Diploma thesis! Thorsten Kemmerer, Michael Lahanas, Dimos Baltas and Nikolaos Zamboglou, DVH computation comparisons using conventional methods and optimized FFT algorithms for brachytherapy, Med. Phys. 27 2343-2356, 2000 Abstract
Last Update 1-3-2004
|
|
from Ns sources, where the i-th source at 
produces a dose 
at the point 
. The total dose distribution the point 
|
Contact - Search - Statistics |
Counter Start 20-4-2004