A feasibility study to estimate optimal rigid‐body registration using combinatorial rigid registration optimization (CORRO)

Abstract Purpose Clinical image pairs provide the most realistic test data for image registration evaluation. However, the optimal registration is unknown. Using combinatorial rigid registration optimization (CORRO) we demonstrate a method to estimate the optimal alignment for rigid‐registration of clinical image pairs. Methods Expert selected landmark pairs were selected for each CT/CBCT image pair for six cases representing head and neck, thoracic, and pelvic anatomic regions. Combination subsets of a k number of landmark pairs (k‐combination set) were generated without repeat to form a large set of k‐combination sets (k‐set) for k = 4,8,12. The rigid transformation between the image pairs was calculated for each k‐combination set. The mean and standard deviation of these transformations were used to derive final registration for each k‐set. Results The standard deviation of registration output decreased as the k‐size increased for all cases. The joint entropy evaluated for each k‐set of each case was smaller than those from two commercially available registration programs indicating a stronger correlation between the image pair after CORRO was used. A joint histogram plot of all three algorithms showed high correlation between them. As further proof of the efficacy of CORRO the joint entropy of each member of 30 000 k‐combination sets in k = 4 were calculated for one of the thoracic cases. The minimum joint entropy was found to exist at the estimated mean of registration indicating CORRO converges to the optimal rigid‐registration results. Conclusions We have developed a methodology called CORRO that allows us to estimate optimal alignment for rigid‐registration of clinical image pairs using a large set landmark point. The results for the rigid‐body registration have been shown to be comparable to results from commercially available algorithms for all six cases. CORRO can serve as an excellent tool that can be used to test and validate rigid registration algorithms.


| INTRODUCTION
Image guided radiation therapy (IGRT) involves imaging patients before and during radiation treatment with the aim of increasing the dose to the tumor while minimizing radiation to healthy tissues. One of the implementations of IGRT requires the ability to register the daily cone beam computed tomography images (CBCT) to the reference planning CT image set. 1 The calculated patient position correction is then used to shift the patient into near perfect alignment with the reference image position prior to beam delivery. Image guided radiation therapy depends extensively on image registration for spatial alignment and quantitative evaluation for changes in the position and size of the target and normal tissues due to weight loss and tumor response. 2 The registration process geometrically aligns two images by finding the minimum of an objective function representing the alignment quality (typically the mean squared error to determine the parameters of a rigid registration transformation matrix). Validation and quantification of image registration quality still poses a very challenging problem in clinical practice due to the lack of an underlying ground truth. Currently, there are four leading methods to assess registration quality: visual inspection, use of fiducials, landmark point sets, and mutual information.
Conventionally, physicians have validated registered images by visually inspecting portal images and diagnostic quality images alongside a planning digital reconstructed radiograph (DRR). 3 This method of assessing image registration quality typically involves the use of an in-field metric (graticule mounted to the MV treatment head) to identify anatomical structures shared between the portal image and DRR, typically the bony anatomy. The accuracy of this method has been reported to be between 5 to 10 mm. 4 This type of registration is subjective and cannot be used for large quantities of data. 5 Fiducial markers on phantoms and patients have also been used to find the ground truth of a rigid image registration. 6 Due to the visibility of gold markers placed in the prostate in kV or MV x-ray imaging, they are used to improve targeting for surgical procedures and radiotherapy treatment to the prostate and to estimate registration error. 7,[8][9][10] Although fixed to the target anatomy, the fiducials are known to drift from the original fixed point either due to anatomy changes over time or detachment of the fiducial. The effect of fiducial relocation may lead to interobserver error associated with the registration. [11][12][13] Also, the number of fiducials that can be fixed at any time is limited. O'Neill et al. 14 described that in a study of 427 patients undergoing intensity modulated radiation therapy using fiducial marker IGRT the intrafraction motion was found to be greater than 2mm for about 66% of their patients.
Expert positioned landmark point pairs have also been used to quantitatively validate registration quality. These landmark points are used to specify the ground truth in correlated images which allows for validation based on the accuracy of the manual or automatic selection of the points. [15][16][17] In this study, a metric for registration quality assessment through the development of a new approach to rigid registration we call combinatorial rigid registration optimization (CORRO). In CORRO, we generate large landmark sets from an expert selected anatomical landmark set using the mathematics of combination without replacement. We present a method for the quantitative measurement of registration quality between CORRO and commercially available rigid registration software based on mutual information and joint entropy thickness for all cases. The machine isocenter is located at the center of the CBCT reconstruction image volume. The planning CT was resampled to the same dimensions in the X/Y plane as the CBCT and the image content was shifted to place the anatomic isocenter at the center of the planning target volume as shown in Fig. 1.

2.C | Rigid-body registration
Given the CBCT (P) and the planning CT (Q), we calculate the transformation Q = T(P) such that the corresponding coordinates in the two images correspond to the same physical location in both images.
gdenote the collection of points in  3 with the same size, with P representing landmark points in the planning CBCT image and Q representing corresponding landmark points in the CT image. The registration problem in three-dimension (3D) consists of finding the transformation that achieves the best match between the corresponding features in P and Q such that the root mean square (RMS) distance d i between corresponding points is minimized. 9,18 The appropriate translation vector is simply the mean displacement between the two sets of points. 8 The aim is to find the errors associated with locating the landmark points. The image registration problem is reduced to a shape analysis problem or to the YORKE ET AL.

| 15
orthogonal Procrustes problem. 8,9,19,20 The Procrustes problem is simply a least square-fitting problem and studies have shown that the calculation of the rotation matrix R is more involved due to the nonlinear condition for a rotation matrix to be orthogonal. If P and Q are replaced with their centroid values then the optimal transformation is represented as This reduces the problem to the orthogonal Procrustes problem where we seek to find the rotation R. The RMS distance to be minimized is termed as the fiducial registration error (FRE). 8,9 Therefore, given four noncoplanar points for a 3D volume the problem of rigid-body registration is to finding a rotation and translation (t) that minimizes the FRE which is represented mathematically as An FRE value of zero means the rigid-body registration is perfect. However, the FRE is expected to me more than zero due to error associated with locating the landmark points.
FRE can be used to characterize the magnitude of the error associated in locating these landmark points. The translation is given by where the bar indicates a mean over i = 1,⋯, k. In this study, no rotation is used though in general it could be used for robotic couch that allows for 6 degree of rotation.

2.D | Generation of large k-Sets
Using the anatomical landmark pairs that have been selected by expert we generate our k-set as follows. A k-combination set was generated as a subset of the landmark set by combining k number of landmark pairs without repeat. The set of all possible unique k-combinations forms a large set of discrete independent trials, which we term the k-set. In this study, three different sizes of the k-combination set were used, k = 4, 8, 12. The mathematics is similar to the problem of determining the number of unique hands of cards from a standard deck of 52.
Let N I f gi ¼ 1,n be the locations for the landmark points selected by the expert. We find the combination without repetition of the elements such that The points are then used to find the k-combination sets for k = 8 and k = 12. CORRO is designed to estimate the optimal registration output for rigid-registration based on a large set of discrete independent trials, which statistically validate our results using the central limit theorem. The central limit theorem of statistics states that if a large sample is drawn from a population, the distribution of the sample mean is approximately normal and the standard deviation decreases as the number of samples increases. 21 The larger the set of combination sets, the smaller the registration error which leads to a better estimate of the optimal registration mean of the rigid registration. This statistical theorem has been tested by Castillo et al.,17 who developed a framework for evaluating deformable images registration.

2.E | Estimating the ground truth
Having generated a large set of k-combination sets for each case we are able to estimate the k-mean of the rigid-registration and its associated error. The k-mean is found by calculating the mean of all the translations given by the k-set registration outputs and the standard deviation (registration error). The results are validated using the central limit theorem evidenced by Fig. 3 for the mean X-translation for sample case 3.
The k-set is a large population of sets of paired points, (kcombination sets) when we draw a large sample from the population of a k-set the distribution of the sample mean approaches a normal distribution, and the standard deviation of the sample mean decreases as the sample size increases seen in Fig. 3(b). The sampling distribution of the sample mean is found by drawing repeated random samples of a given size n from the population of k-set for a particular case which has a population mean and a population standard deviation. The mean of all the random samples drawn from the population approximate the population's mean.
The Gaussian distribution for one case for k-set of k = 4 is shown in Fig. 3(a) and the standard deviation, which is proportional to 1= ffiffi ffi n p , where n is the sample size, is demonstrated in Fig. 3(b). It should be noted that for this particular case the estimated registration (population mean) for about 398 000 k-set is 5.6816 AE 0.0771 mm for x-translations and the mean of the sample mean is also 5.7048 AE 0.0718 mm as expected by the central limit theorem.

2.F | Joint histogram and joint entropy
Information theory is the theoretical foundation on which information content can be quantified. 22 Entropy measures the randomness or the content of information. Entropy is related to the definition of entropy in thermodynamics, which measures the molecular disorder, or the randomness of a system.
Let X be a discrete random variable over χ and probability density function f (x), x ∈ χ. The entropy of X is defined as where f(x)∈ [0,1], ∑ | 17 A probability of zero does not contribute to the entropy and as a measure of the average uncertainty in X, the entropy is always nonnegative and indicates the number of bits on average required to describe the random variable. The higher the entropy the more information the variable contains. Consequently, the concept of entropy can be extended to two or more variables. The joint entropy of a pair of random variables H (X, Y) of the planning CT (x) and the registered CBCT (y) with joint probability distribution is defined as.
H (X, Y) is always at least equal to the entropies of X and Y alone.
That is, adding a new variable cannot reduce the existing uncertainty.
The idea of calculating joint entropy for a pair of images is to find the number of possible grey scale values in each image using the 2D joint histogram f (x, y). If the images are perfectly aligned then the histogram is extremely focused. Misalignment between the images means that there is more dispersion meaning there is higher entropy, as the entropy is the measure of the randomness. In effect images are registered when the coordinate in one image space relates that of the other image to minimize the joint entropy and the randomness is minimized.
The opposite is true for unregistered images.

3.A | Image Data and landmark pairs
Using Eq. (6) the joint entropy is calculated and a summary of the joint entropy for all six cases is presented in  Table 2 shows the number of landmarks pairs that were selected for each CT-CBCT case. It should be noted that not every discrete independent combination would give an acceptable k-set. In this result the same k-set was reported for k = 4, 8, and 12. Figure 5 shows the combined joint histogram distribution from CORRO (red), AdaPT Insight (Green) and MIM (Blue) for case 3.

3.B | Trials of Landmark Sets and Rigid-Body Registration
Large k-sets for k = 4, 8 & 12 were generated and used to solve for the rigid-body registrations the optimal registration was found by calculating the population mean from the k-set for a k = 4,8 &12; the standard deviation was found to be the registration error. Table 1 shows the total number of independent trials that were successfully used to perform the registration for each of the cases. Applying the mean registration for combination k = 4 for case 2 to the CBCT the output registration is demonstrated by fusing the planning CT and CBCT to see how well the rigid registration performed. Figure 4 gives a pictorial view of the registration result displayed in a checkerboard form and compared to the registration from the two commercially available registration software packages AdaPT Insight and MIM respectively. Landmark features are well-aligned in CORRO compared to commercial software.

3.C | Validating Joint Entropy Results
Computing the joint entropy of the planning CT and CBCT combination by targeting only the tumor regions for all six cases also This can be seen in Table 1 for case 3 where the joint entropy for   CORRO and AdaPT Insight are much closer in value than CORRO and MIM or AdaPT Insight and MIM.
The joint entropy for all six cases is tabulated in Table 2 where the results from CORRO are compared to those from AdaPT Insight and MIM software. For all six cases the joint entropy from CORRO were very comparable with very slight differences. Fig. 6 also shows the joint entropy for the 30,000 translations from case 5 for k = 4.
The joint entropy for each output was calculated and plotted against the root mean squared distance. The joint entropy is calculated for the image using each member in the k-set. For example, if there are 190,000 members in the k-set there will be a total of 190 000 joint entropies calculated. Figure 6 is calculated using 30 000 k- combination sets the mean translations in the X, Y, Z plane are 9.8, 0.76, and 0.96-mm respectively. When joint entropy was plotted against the X translation the minimum joint entropy was found to exist at 9.83 mm in the X direction. The minimum joint entropy was found to exist at 0.8 mm pixels when plotted against the Y translation and was found to exist at 1.0 mm when plotted against the Z translation. Plotting the joint entropy vs the RMS value of X, Y, Z translation as seen in Fig. 6(d) the minimum joint entropy was found to exist exactly at the root mean square distances of X, Y, Z i.e. ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi X 2 mean þ Y 2 mean þ Z 2 mean r = 9.88 mm which is expected.

| DISCUSSION
In order to objectively evaluate clinical results of an image registration algorithm, better attention needs to be given to the selection of appropriate reference points, with which an algorithm performance can be compared. In their work, Wong et al. 24  The registered images are validated using the root mean squared (RMS) difference of the intensities of the two images, median-absolute deviation of the intensity difference, and maximum intensity differences. 25 The results were validated using the central limit theorem and found that, the estimated error associated with rigid registration reduces as we increase the sample size for a particular k-set as seen in Fig. 3. The results were also validated using joint entropy calculations. A smaller joint entropy value shows there is less variation in the pixel intensities and the similarity between the images is high implying a better registration. Joint Entropy is a common method employed to validate the quality of radiographic images, analyze the similarity of images, and also validate registration output between multimodal image registrations. 26,27 The registration results from CORRO using the mean as the optimal registration were compared to commercially available registration software using joint histogram and joint entropy analysis. A comparison between the joint entropy values of the rigid registration and those from the commercially available software were very closely related as seen in Table 2. The results of k = 4 gave lower joint entropy values when compared to k = 8 and 12 and the commercial software indicating slightly better registration results using k = 4.
This was not surprising because k = 8 and 12 are more likely to have mismatched points that would affect the accuracy of the registration. 28 In 3D any four points can be mapped to another four points. 29 However, more landmarks result in transformations with nonuniform biased sampling of the image. Coste 30 shows that using many landmarks could result in odd transformations of the image grid. Our data comparing CORRO for k = 4 and k = 8 and k = 12 seem to support such a hypothesis k = 4 demonstrates superior registration results in comparison to k = 8 and 12. Hence k = 4 can be best used to determine the optimal registration results.
Lastly, the joint entropy for all the translations plotted against the root mean square (RMS) distances of the output translations (X, Y, Z). The root mean square distance is given by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Xmean 2 þ Ymean 2 þ Zmean 2 p . We found the minimum joint entropy value to be exactly equal to the RMS value. Signifying that the mean translation in X, Y, Z are equivalent to the optimal rigid registration translation outputs as determined by the calculated minimum joint entropy. It should be noted that the landmark selection process is time consuming. Hence this study will not be suitable for clinical routine. This work was carried out to demonstrate the feasibility of finding the optimal rigid transformation using large sets of landmark points; in the future automated techniques could be used to identify the landmark points.

| CONCLUSION
The results of this study show that we have developed a methodology called CORRO that allows us to estimate optimal alignment F I G . 5. Combination of Joint Histogram distributions from combinatorial rigid registration optimization (Red), AdaPT Insight (Green) and MIM (Blue) for case 3.
for rigid-registration of clinical image pairs using a large set of discrete independent trials of large landmark points combination sets.
The viability of this method has been demonstrated with six cases: two pelvic, two head and neck, and two thoracic. The results for the rigid-body registration have been shown to be comparable to results from commercially available algorithms for all six cases.
CORRO can serve as an excellent metric for registration quality assessment. The next step in this project is using CORRO to create a statistically characterized reference data set for 58 pelvic cases that are being made available at the cancer imaging archives (https://wiki.cancerimagingarchive.net/display/Public/Pelvic+Refere nce+Data).

ACKNOWLEDGMENTS
We thank Beaumont Hospital Royal Oak MI, for their support. We also thank C. K. Yorke, Ph.D., B. J. Roth Ph.D., and G. C. McDonald Ph.D., D.Sc. for their assistance and enlightening comments.

CONFLI CT OF INTEREST
Authors declare no financial relationships that may lead to a conflict of interest.