**Fast calculation of the exact radiological path for a 3D CT array**

*Robert L. Siddon*

**ABSTRACT**

Ready availability has prompted the use of computed tomography (CT) data in various applications in radiation therapy. For example, some radiation treatment planning systems now utilize CT data in heterogeneous dose calculation algorithms. In radiotherapy imaging applications, CT data are projected onto specified planes, thus producing *radiographs*, which are compared with simulatorradiographs to assist in proper patient positioning and delineation of target volumes. All these applications share the common geometric problem of evaluating the radiological path through the CT array. Due to the complexity of the three-dimensional geometry and the enormous amount of CT data, the exact evaluation of the radiological path has proven to be a time consuming and difficult problem. This paper identifies the inefficient aspect of the traditional exact evaluation of the radiological path as that of treating the CT data as individual voxels. Rather than individual voxels, a new exact algorithm is presented that considers the CT data as consisting of the intersection volumes of three orthogonal sets of equally spaced, parallelplanes. For a three-dimensional CT array of N^{3} voxels, the new exact algorithm scales with 3N, the number of planes, rather than N^{3}, the number of voxels. Coded in FORTRAN-77 on a VAX 11/780 with a floating point option, the algorithm requires approximately 5 ms to calculate an average radiological path in a 100^{3} voxel array.

**INTRODUCTION**

In radiation therapy applications, computer tomography (CT) data are utilized in various dose calculation and imaging algorithms. For example, some radiation treatment planning systems now utilize two-dimensional CT data for pixel-based heterogeneous dose calculations. Other systems forward project three-dimensional CT data onto specified planes, thus forming *radiographs*, which are compared with simulator radiographs to assist in proper patient positioning and delineation of target volumes. All such applications, whether in inhomogeneity calculations or imaging applications, essentially reduce to the same geometric problem: that of calculating the radiological path for a specified ray through the CT array.

Although very simple in principle, elaborate computer algorithms and a significant amount of computer time is required to evaluate the exact radiological path. The amount of detail involved was recently emphasized by Harauz and Ottensmeyer, who stated that even for the two-dimensional case, their algorithm for calculating the exact radiological path grew more and more unwiedly and time consuming, while remaining unreliable. For three dimensions, they concluded that determining the exact radiological path is not viable. This paper describes an *exact*, efficient, and reliable algorithm for calculating the radiological path through a three-dimensional CT array.

Denoting a particular voxel density as ρ(i,j,k) and the length contained by that voxel as *l*(i,j,k), the radiological path may be written as

(1) *d* = Σ_{i}Σ_{j}Σ_{k}*l*(i,j,k)ρ(i,j,k)

Direct evaluation of Eq.(1) entails an algorithm which scales with the number of terms in the sums, that is, the number of voxels in the CT array. The following describes an algorithm that scales with the sum of linear dimensions of the CT array.

**METHOD**

Rather than independent elements, the voxels are considered as the intersection volumes of orthogonal sets of equally spaced, parallel planes. Without loss of generality, Fig. 1 illustrates the two-dimensional case, where pixels are considered as the intersection areas of orthogonal sets of equally spaced, parallel lines. The intersections of the ray with the lines are calculated, rather than intersections of the ray with the individual pixels.

Determining the intersections of the ray with the equally spaced, parallel lines is a particularly simple problem. As the lines are equally spaced, it is only necessary to determine the first intersection and generate all the others by recursion. As shown on the right illustration of Fig. 1, the intersections consist of two sets, one set for the intersections with the horizontal lines (closed circles) and one set for the intersections with the vertical lines (open circles). Comparing the left and right illustrations of Fig. 1, it is clear that the intersections of the ray with the pixels is a subset of the intersections with the lines. Identifying that subset allows the radiological path to be determined. The extension to the three-dimensional CT array is straightforward.

The ray from point 1 to point 2 may be represented parametrically as

(2) X(α) = X_{1}+α(X_{2}-X_{1}), Y(α) = Y_{1}+α(Y_{2}-Y_{1}), Z(α) = Z_{1}+α(Z_{2}-Z_{1})

where the parameter α is zero at point 1 and unity at point 2. The intersections of the ray with the sides of the CT array are shown in Fig. 2.

If both points 1 and 2 are outside the array [Fig. 2(a)], then the parametric values corresponding to the two intersection points of the ray with the sides are given by α_{min} and α_{max}. All intersections of the ray with individual planes must have parametric values which lie in the range (α_{min},α_{max}). For the case illustrated in Fig 2(b), where point 1 is inside the array, the value of α_{min} is zero. Likewise, for Fig. 2(c), if point 2 is inside, then α_{max} is one. For both points 1 and 2 inside the array [Fig. 2(d)], then α_{min} is zero and α_{max} is one. The solution to the intersection of the ray with the CT voxels follows immediately: Determine the parametric intersection values, in the range (α_{min},α_{max}), of the ray with each orthogonal set of equally spaced, parallel planes. Merge the three sets of parametric values into one set; for example, merging of the sets (1,4,7), (2,5,8), and (3,6,9) results in the merged set (1,2,3,4,5,6,7,8). The length of the ray contained by a particular voxel, in units of the ray length, is simply the difference between two adjacent parametric values in the merged set. For each voxel intersection length, the corresponding voxel indices are obtained, and the products of the length and density are summed over all intersections to yield the radiological path. A more detailed description of the algorithm is given in the following section.

**ALGORITHM**

For a CT array of (N_{x}-1, N_{y}-1, N_{z}-1) voxels, the orthogonal sets of equally spaced, parallel planes may be written as

(3) X_{p}(i) = X_{p}(0) + id_{x}, Y_{p}(j) = Y_{p}(0) + jd_{y}, Z_{p}(k) = Z_{p}(0) + kd_{z},

where i, j, and k are integers and d_{x}, d_{y}, and d_{z} are the distances between the x, y, and z planes, respectively. The quantities d_{x}, d_{y}, and d_{z} are also the lengths of the sides of the voxel. The parametric values α_{min} and α_{max} are obtained by intersecting the ray with the sides of the CT array. As shown in Fig. 2(d) both 1 and 2 are assumed to be located on one of the three equally spaced, parallel planes. This is also assumed to be the case if 1 and/or 2 are located outside the CT array. From Eqs. (2) and (3), the parametric values corresponding to the sides are given by the following:

(4) If(X_{2}-X_{1})≠0, α_{x}(0)=(X_{p}(0)-X_{1})/(X_{2}-X_{1}), α_{x}(N_{x}-1)=(X_{p}(N_{x}-1)-X_{1})/(X_{2}-X_{1})

with similar expressions for α_{y}(0), α_{y}(N_{y}-1), α_{z}(0), and α_{z}(N_{z}-1). If the denominator (X_{2}-X_{1}) in Eq. (4) is equal to zero, then the ray is perpendicular to the x axis, and the corresponding values of α_{x} are undefined, and similarly for α_{y} and α_{z}. If the α_{x}, α_{y}, or α_{z} values are undefined, then those values are simply excluded in all the following discussion.

In terms of the parametric values given above, the quantities α_{min} and α_{max} are given by

(5.1) α_{min}=max{0, min[α_{x}(0),α_{x}(N_{x}-1)], min[α_{y}(0),α_{y}(N_{y}-1)], min[α_{z}(0),α_{z}(Nz-1)]},

(5.2) α_{max}=min{0, max[α_{x}(0),α_{x}(N_{x}-1)], max[α_{y}(0),α_{y}(N_{y}-1)], max[α_{z}(0),α_{z}(Nz-1)]},

where the functions *min* and *max* select from their argument list, the minimum and maximum terms, respectively. If α_{max} is less than or equal to α_{min}, then the ray does not intersect the CT array.

From all intersected planes, there are only certain intersected planes which will have parametric values in the range (α_{min},α_{max}). From Eqs. (2), (3), and (5), the range of indices (i_{min},i_{max}), (j_{min},j_{max}), and (k_{min},k_{max}), corresponding to these particular planes, are given by the following: