Real-time estimation of prostate tumor rotation and translation with a kV imaging system based on an iterative closest point algorithm

There are many ways to represent rotation, such as: Euler angles, Gibbs vector, Cayley–Klein parameters, Pauli spin matrices, axis and angle, orthonormal matrices and Hamilton's quaternions (Korn and Korn 1968, Stuelpnagle 1964). In this study, Euler angles were used as it is more convenient for real-time tracking on a linac system; rotational variations about the RL, SI and AP axes are the Euler angles θ_RL, θ_SI and θ_AP, respectively.

The estimation of rotation matrix and translation vector parameters (also called pose estimation) is a common problem in computer vision. One of the first studies was by Thompson (1958) who developed a method that calculates rotation when three data points are available. The method of Thompson (1958) depends on selective neglect of the extra constraints when all coordinates of three points are known. Schut (1960) proposed an alternative method based on unit quaternion and a set of linear equations. However, neither method handles more than three data points nor they use all the information available from the three points. More recently, Arun et al (1987) developed a numerical least squares fitting method using only two data points. Their algorithm was based on the singular value decomposition (SVD) of a 3 × 3 matrix. A closed form solution to the least squares problem for three or more points was developed by Horn (1987). Deriving the solution was simplified by using unit quaternion to represent rotation. However, due to reflection, the Arun and Horn methods fail to determine the correct rotation matrix when the data is extremely noisy. These problems with reflection were addressed by Umeyama (1991) by changing the matrix of the solution.

In this study, we implemented the Umeyama method and the ICP algorithm (figure 1) to estimate rotation and also examine the maximum error level according to this estimation (Chen and Medioni 1991). Figure 1 shows the steps for the ICP algorithm in more detail.

Zoom In Zoom Out Reset image size

A total of 11 748 images acquired during one fraction of the treatment of ten prostate cancer patients with radiopaque prostate markers were used in this study (Ng et al 2012). Images during VMAT treatment were acquired by a kV x-ray imager mounted perpendicularly to the radiation treatment beam source (OBI, Varian Medical Systems). All treatments were dual-arc VMAT with similar treatment times and the number of images per patient was similar. The exposure parameters used were 125 kVp, 80 mA and 13 ms (a standard pelvic cone beam computed tomography scan setting) with a 6 × 6 cm² field size (no filter was used). The imaging frequency was 5 or 10 Hz, depending on the desired image quality. During treatment, to reduce scatter from the MV source, the kV source to detector distance was set to 180 cm (compared with 150 cm for CBCT). The estimated imaging dose for kilovoltage intrafraction imaging is 185 mGy over the course of a conventionally fractionated prostate IMRT treatment. More details of the dose estimates are given in Crocker et al (2012).

As the gantry rotates around the patient during treatment, two dimensional projections of the prostate are acquired at time t_k by a kV imager. The fiducial markers are segmented and extracted using in-house software (Fledelius et al 2011b). Figure 2 is a screenshot showing a kV image with three radiopaque markers and the segmented location of the markers as determined by our in-house software.

Zoom In Zoom Out Reset image size

To determine the three dimensional location of the markers from a series of two dimensional segmented markers, a three dimensional Gaussian probability distribution function (PDF) is assumed for each marker (Poulsen et al 2008b). The influence of the PDF does not affect the comparison because in both methods (three dimensional and ICP algorithm) the calculation of the PDF is the same. However, the use of the PDF itself (and errors in marker segmentation) may induce uncertainty in the rotation estimation. The covariance matrix and mean marker location used by the PDF are determined by fitting the PDF to a series of two dimensional marker locations by the maximum likelihood estimation. The three dimensional locations of the markers are determined from the three dimensional PDF using the method of Poulsen et al (2008a).

Let us assume that we have two sets of data points for the fiducial markers, a source point at time t₁ and a target point at time t_k. The source points are determined during a pretreatment arc that is necessary to build the PDF needed for the translation (and rotation) estimation prior to turning the treatment beam on (Ng et al 2012). The target point at time t_k is obtained by rotating the source point, t₁, around an axis and/or translating, the source point t₁, along an axis with an unknown transform function. As implemented the rotation matrix of the ICP algorithm is computed about the center of mass of the markers after the shift of the translation T. The set of source points, one for each marker, are denoted by Q = (q₁, q₂, ..., q_N) and the set of target points, one for each marker, are denoted by P = (p₁, p₂, ..., p_N), where N is the number of fiducial markers. We assume that all markers are translated and rotated by the same amount, so that

$$\begin{equation} Q( P ) = sRP + T + n \end{equation} \tag{ 1 }$$

where R is the rotation matrix, T is the translation vector, s is the scaling factor and n is the noise (assuming the noise is Gaussian). We separate the rotation matrix into its components about each axis:

$$\begin{equation} R = R_{{\rm RL}} ( {\theta _{{\rm RL}} })R_{{\rm SI}} ( {\theta _{{\rm SI}} } )R_{{\rm AP}} ( {\theta _{{\rm AP}} }) \end{equation} \tag{ 2 }$$

$$\begin{equation} R_{{\rm RL}} ( {\theta _{{\rm RL}} } ) = \left[ {\begin{array}{@{}ccc@{}} 1 & 0 & 0 \\ 0 & {\cos ( {\theta _{{\rm RL}} })} & { {-} \sin ( {\theta _{{\rm RL}} } )} \\ 0 & {\sin ( {\theta _{{\rm RL}} })} & {{\rm cos}( {\theta _{{\rm RL}} })} \\ \end{array}} \right] \end{equation} \tag{ 3 }$$

$$\begin{equation} R_{{\rm SI}} ({\theta_{{\rm SI}}}) = \left[ {\begin{array}{@{}ccc@{}} {\cos ({\theta _{{\rm SI}} })} & 0 & {\sin ( {\theta _{{\rm SI}} })} \\ 0 & 1 & 0 \\ { {-} \sin ( {\theta _{{\rm SI}} } )} & 0 & {\cos ( {\theta _{{\rm SI}}})} \\ \end{array}} \right] \end{equation} \tag{ 4 }$$

$$\begin{equation} R_{{\rm AP}} ( {\theta _{{\rm AP}} }) = \left[ {\begin{array}{@{}ccc@{}} {{\rm cos}( {\theta _{{\rm AP}} } )} & { - {\rm sin}( {\theta _{{\rm AP}} } )} & 0 \\ {{\rm sin}( {\theta _{{\rm AP}} } )} & {{\rm cos} ( {{\rm }\theta _{{\rm AP}} {\rm }} )} & 0 \\ 0 & 0 & 1 \\ \end{array}} \right]. \end{equation} \tag{ 5 }$$

The ICP algorithm iteratively performs matching for every data point based on the nearest neighbor algorithm (Cover and Hart 1967), see step 1 in figure 1. In our context, this step is not normally needed because it is easy to identify which marker the data point belongs to. The nearest neighbor analysis is very important when two or more of the fiducial markers are close to each other and could be mistakenly identified. The three fiducial markers in most of the data studied in this research were placed sufficiently far apart to avoid mismatching. In clinical practice, either sufficient marker separation would be required or robust methods to account for the issues with overlapping markers would need to be developed.

For step 2 in figure 1, we plan to estimate the rotation and transformation matrices, R and T. To do this, we minimize the sum of the squared distances of target to source points (Besl and McKay 1992). This forms the error function f(R, T) ≈ n (n is the noise in equation (1)) that we wish to minimize and is mathematically expressed as

$$\begin{equation} \mathop {\rm arg\min }\limits_{R,T} f ( {R,T,Q,P} ) = \frac{1}{N} \sum _{i = 1}^N ( {q_i - ( {Rp_i + T} )} )^2 . \end{equation} \tag{ 6 }$$

We estimate R and T using SVD (Umeyama 1991) which is described in more detail in the appendix. Note that in equation (6) we have set the scale factor, s, to 1. The scale factor will be calculated based on the area of triangles after calculating rotation and translation. The scale factor is essentially the relative area of the triangle made by the three markers. A variation in the scale could indicate errors in marker segmentation or indicate prostate deformation. During an individual fraction marker segmentation errors are more likely than prostate deformation, so we can use the scale factor as a threshold to remove suspected marker segmentation errors.

Using our estimates of R and T, we iteratively update the values using the following algorithm (step 3–5 in figure 1). We use two intermediate matrices, R_it which we initially set to the 3×3 identity matrix, and T_it which we initially set to the 3×1 zero vector. These two matrices are updated according to

$$\begin{equation} R_{it} ({k + 1}) = RR_{it} (k) \end{equation} \tag{ 7 }$$

$$\begin{equation} T_{it} ( {k + 1} ) = RT_{it} ( K )P_k + T \end{equation} \tag{ 8 }$$

where k is the iteration number; in this case we terminate after a maximum of five iterations. In step 4 of figure 1, the new point P_{k + 1} is calculated from the new update of the following equation:

$$\begin{equation} P_{k + 1} = R_{it} ( {k + 1} )P_k + {\rm repmat} ( {T_{it} ( {k + 1} )1,3} ) \end{equation} \tag{ 9 }$$

where repmat (T_it(k + 1)1, 3) is a matrix consisting of a 3 × 1 tiling of copies of T_it(k + 1). Now that we have a new value for P_{k + 1}, we calculate new estimates for R and T using equation (6), i.e. the new estimates, R_{k + 1}, and T_{k + 1} are found with

$$\begin{equation*} arg\mathop {\min }\limits_{R_{k + 1} = R,T_{k + 1} = T} f ( {R,T,Q,P_{k + 1} }). \end{equation*}$$

Convergence usually occurs in fewer than 5 iterations (see figure 3 for an example of the error for 21 iterations). In our implementation we perform at most five iterations and the final value of P_{k + 1} is used to estimate the root mean square error (RMSE) of the matrix transformation using equation (10)

$$\begin{equation} {\rm RMSE} = \sqrt {\sum \limits_{i = 1}^N \frac{{( {p_{k + 1} ( i ) - q( i )} )^2 }}{N}} . \end{equation} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

The final step is to calculate the Euler angles, θ_RL, θ_SI, θ_AP, from the rotation matrix, R, using a nonlinear least square fitting algorithm (Klumpp 1976, Weisstein 2013) for the new fiducial marker positions.

A program has been written in C# to show the 3D trajectory of the markers as calculated in Poulsen et al (2008b) with the addition of the rotation calculations as presented in this study. Figure 4 shows the main screen of the software which shows, in real-time, the latest kilovoltage image with the segmented 2D marker locations, a graphical view of the markers in the middle, the 3D trajectory of the markers and the rotation of the markers.

Zoom In Zoom Out Reset image size

In order to evaluate the method, we simulated the rotation of three fiducial markers based on assuming the 3D positions were perfectly estimated. In each simulated image, the triangle formed by the three fiducial markers was rotated in each subsequent image by θ_RL = 0.1°, θ_SI = 0.2°, θ_AP = 0.3°. Figure 5 shows the value of the rotation around RL, rotation around SI and rotation around AP sequentially.

Zoom In Zoom Out Reset image size

The mean of the RMSE between the source points and transformed target points (equation (10)) were calculated up to 30°. Figure 5 shows the linear incremental rotation of up to 30° for three directions. The maximum error up to 30° was 0.007° with a mean value of 0.004°. Rotations of 30° are above the rotations that we observed in the next section for prostate cancer patients.

In this section, we extended our analysis to prostate cancer patients where we calculated the real-time prostate rotation and translation of the prostate for ten prostate cancer patients. In figures 6(a) and (b) an example of the time series of rotation and translation is illustrated for a prostate cancer patient. Rotations of over 7° are observed for this patient. The translation and rotation data are relatively noisy because of the noisy data in the 2D segmented location of the markers. However, a clear trend in the rotation data is visible in figure 6(a).

Zoom In Zoom Out Reset image size

Figure 6(c) shows variations in the scale factor, s, for the same patient. The scale factor was calculated based on the changes in the area of the triangles formed by the three fiducial markers relative to the area of the triangle in the first image. Figure 7 gives a graphical view of the triangle formed by the markers. Figure 8 shows the variation of the scale factor for the whole dataset (11 748 images). The correlation coefficient of the template based marker segmentation (Fledelius et al 2011b) and scale factors were implemented to detect and remove noise in the rotation calculations. In this study, the threshold for correlation was 0.65 and the threshold for the scale factor set for values with more than a 5% change in size (1 ± 0.05). Less than 2% of the data were removed based on these thresholds. Figure 9 shows the histogram of the translation and rotation parameters. For rotation, the center of θ_RL was around 1°, while θ_SI and θ_AP were centered around 0°. For translation, the SI and RL histograms are centered close to 0 mm, while the center of the AP histogram is around 0.3–0.4 mm. Figure 10 shows the pair-to-pair datapoints and the correlation of the variables to each other in the upper right hand triangle. The maximum correlation was between SI and AP translations with a correlation coefficient of 0.66. This is consistent with the results from Poulsen et al (2008b). The correlations between θ_RL, θ_SI and θ_AP also had values higher than 0.1 with a negative correlation between θ_RL and θ_AP. Figure 10 also shows that the correlation between θ_SI and θ_AP was almost two times larger than correlation between θ_RL and θ_SI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In table 1 the value of SD of rotation and translation is given in millimeters and degrees, respectively. The maximum value of SD is θ_RL with 2.3° of rotation.

Table 1. Standard deviation of rotation and translation.

$\theta _{{\rm RL}} ^\circ$	$\theta _{{\rm SI}} ^\circ$	$\theta _{{\rm AP}} ^\circ$	RL (mm)	SI (mm)	AP (mm)
2.26	0.89	0.72	0.32	0.53	1.13

Figure 11 shows the RMSE values between the source points and transformed target points (equation (10)) of the rotation with the ICP and translation only. The x-axes show the distance error of the transformed points using the new rotation matrix R and translation vector T. The mean, median and SD of the errors is given in table 2. The ICP algorithm had a sub-0.2 mm mean, median and SD. These values were significantly smaller compared to those of the translation only algorithm.

Zoom In Zoom Out Reset image size