Automatic 3D MRI-Ultrasound Registration for Image Guided Arthroscopy

Abstract

Registration of partial view intra-operative ultrasound (US) to pre-operative MRI is an essential step in image-guided minimally invasive surgery. In this paper, we present an automatic, landmark-free 3D multimodal registration of pre-operative MRI to 4D US (high-refresh-rate 3D-US) for enabling guidance in knee arthroscopy. We focus on the problem of initializing registration in the case of partial views. The proposed method utilizes a pre-initialization step of using the automatically segmented structures from both modalities to achieve a global geometric initialization. This is followed by computing distance maps of the procured segmentations for registration in the distance space. Following that, the final local refinement between the MRI-US volumes is achieved using the LC${}^2$ (Linear correlation of linear combination) metric. The method is evaluated on 11 cases spanning six subjects, with four levels of knee flexion. A best-case error of 1.41 mm and 2.34${}^{\circ }$ and an average registration error of 3.45 mm and 7.76${}^{\circ }$ is achieved in translation and rotation, respectively. An inter-observer variability study is performed, and a mean difference of 4.41 mm and 7.77${}^{\circ }$ is reported. The errors obtained through the developed registration algorithm and inter-observer difference values are found to be comparable. We have shown that the proposed algorithm is simple, robust and allows for the automatic global registration of 3D US and MRI that can enable US based image guidance in minimally invasive procedures.

1. Introduction

Knee arthroscopy is a minimally invasive surgery (MIS) commonly used to treat intra-articular knee disorders [1]. The procedure is performed by inserting an arthroscope and a surgical tool into the knee joint through small incisions. The setup allows the surgeon to inspect the joint through the arthroscopic view displayed on a screen while altering the leg position to allow the scope to capture the required tissues. Despite offering several clinical advantages, arthroscopic procedures face long-standing challenges such as physically-demanding ergonomics during patient manipulation, lack of depth perception, limited field of view, and counter-intuitive hand-eye coordination between scope and surgical instruments.Even with skilled operator manipulation, a complex procedure like this can cause unintended trauma and post-operative complications [2,3,4].

To overcome these setbacks, an image guidance platform which uses both intra-operative and pre-operative imaging, in particular US and MRI, is proposed in [5]. This image guidance platform enables the surgeon to track and visualize both the target tissue distribution and the position of the surgical tool(s) in real-time. Intra-operative US is currently the only non-invasive, non-ionizing, imaging modality compatible with the knee arthroscopy surgical setup that can offer real-time imaging for guidance during the surgery [6,7]. However, US can only offer a limited field of view [8], with complex image interpretation. On the other hand, MRI provides a superior soft-tissue contrast, high resolution, as well as field of view covering the entire surgical site (total knee). Pre-operative MRI allows the surgeon to derive and use the information missing from the intra-operative US, enabling a comprehensive view of the bony and soft tissues. Therefore, the combination of both the modalities would synergistically integrate high contrast, large field-of-view, semantically labelable pre-operative MRI data into the fast, intra-operative, partial view 3D US and enable image-guided robotic surgery or MIS. This fusion not only provides a better visualization of structures of interest but also safely enhances the level of automation of MIS. This requires a quick and automatic registration of intra-operative US with pre-operative MRI.

Though relatively new for arthroscopy, multi-modal registration of US and MRI has been previously studied in the field of brain surgery where intra-operative US is used to track brain shift and tumor resection. The MICCAI Challenge CuRIOUS provided a public platform to benchmark MRI-US registration on the RESECT dataset [9]. A total of six fully automated registration algorithms are presented in this challenge [10]. These algorithms were tested on 12 patient volumes and the least registration error, which is 1.57 ± 0.96 mm, was reported by the team ImFusion who performed rigid registration using the LC${}^2$ metric [11]. Another recent work [12] proposed combined segmentation based MRI–US registration approach in image-guided neurosurgery. Completely automatic registration of US and MRI volumes is a complex problem with no satisfactory solution so far [11]. It is primarily because of the entirely different principles used to capture the information in the two modalities [11]. Hence, the direct registration of the two modalities is not easily achievable in the image intensity domain using traditional metrics like the sum of squared distances, mutual information [13] or correlation [14]. Gradient-based methods have also been used, but they tend to discard valuable MRI intensity information [15]. Specifically targeting the multi-modal registration, a similarity measure, LC${}^2$, was proposed in [11], which is invariant to local changes with respect to brightness, contrast, etc. The metric was initially designed for US-CT [16] and later adapted to US-MRI registration.

The other challenging aspect of automating the MRI-US registration is the initialization. In standard clinical practice, close initialization is obtained by selecting 3D landmarks. Defining these landmarks requires not only tedious manual intervention but also keen knowledge of the targeted anatomy and clear geometrical understanding, especially in the case of 3D US, where there is limited volumetric coverage. Hence, to automate the registration process, instead of manually identifying landmarks in 3D, we show that obtaining coarse segmentations and using them for global initialization is a much more convenient alternative [17].

Segmentation-based initialization has been studied for prostate biopsy in [18], where registration was done between trans-rectal US and the corresponding MRI volumes. Even so, the target organ complexity or limited field-of-view of US could demand further regularization of initialization. Hence, we introduce a strong two-level initialization procedure which utilizes structural information of the tissues segmented followed by the distance maps, which were proven to be useful for correspondence estimation [18,19,20].

In this work, we present the first attempt to automatically register intra-operative 4D US to pre-operative 3D MRI as a first step towards the feasibility of autonomous or robotic-assisted knee arthroscopy. Initialization of registration requires strong anatomical priors, which we draw from automatic segmentation of structures from both the modalities. Given the cross-modal nature and partial view obtained in the intra-operative modality (as compared to the view available in pre-operative MRI), an alternative domain that utilizes both the computed priors is adopted—the Euclidean distance transform of segmented masks. This choice addresses the dual needs of working with cross-modal registration, and partial views. The registration of segmented masks in distance space gives us the initial transform for the registration of original volumes. A globally convergent similarity metric, LC${}^2$ [11], is used for this registration.

2. Data

2.1. Ultrasound Data

Ultrasound volumes were acquired by an experienced orthopaedic surgeon F.S. using a Philips EpiQ7 US system and a Philips VL13-5 US probe. SonoCT real-time compound imaging technology and XRES image processing were selected on the US system to enhance image quality. Image acquisition settings have been optimized by a US specialist to work best for knee structure visualization: 13 MHz probe frequency, 4 cm penetration depth and far field focus. Given the restriction of movement due to the presence of surgical tools through the medial/lateral patellar portals, the US probe is positioned on the patellar tendon and parallel to the principal axis of tibia, as shown in Figure 1. 4D US sequences were acquired simulating several possible surgical scenarios for the four flexion angles, 0${}^{\circ }$, 30${}^{\circ }$, 60${}^{\circ }$ and 90${}^{\circ }$ [21]. In this paper, we use 4D US static sequences which were acquired when the leg is fixed at one of the reference knee flexion angles and the probe positioned on the patellar tendon, with the scanning plane parallel to the tibial shaft.

The neutral leg position, where the femoral shaft and tibial tubercle were aligned, was designated as 0${}^{\circ }$ flexion. The flexion angles were then achieved by bending the leg such that the angle between femur and tibia with the neutral position was half the total angle. For angles other than 0${}^{\circ }$, the US scans were recorded with the leg placed on the custom made supporting cushion. Eleven 4D US static sequences comprising 61 3D US volumes of the anterior knee portion were collected from six healthy subjects at different flexion angles.

Table 1. Scans acquired at different angles of flexion for different subjects.

0${}^{\circ }$	30${}^{\circ }$	60${}^{\circ }$	90${}^{\circ }$
Volunteer 1	Volunteer 3	Volunteer 3	Volunteer 3
Volunteer 2	Volunteer 6	Volunteer 6	Volunteer 6
Volunteer 4
Volunteer 5
Volunteer 6

shows the subject wise angle information. Each 3D volume is of 272 × 510 × 256 voxels, corresponding approximately to a size of 4 cm × 4 cm × 5 cm, i.e., each voxel has the dimensions of about 0.14 mm × 0.09 mm × 0.21 mm. The Queensland University of Technology Ethics Committee (No. 1700001110) granted the ethics approval for the US and MRI data collection. The study objective and measurement procedure were explained to each volunteer and their written informed consent is procured.

2.2. MRI Data

MRI images were collected from the same volunteers from whom the US volumes were acquired. The volunteers’ knees were imaged using a 3T MRI system (Siemens Magnetom 3T Prisma, Erlangen, Germany). The MRI scans were acquired with knee positioned at 0${}^{\circ }$, 30${}^{\circ }$, 60${}^{\circ }$ and 90${}^{\circ }$ knee flexion for each volunteer. For the 0${}^{\circ }$ knee flexion, the scan is acquired with the volunteer in the supine position. For the other knee flexion angles, the volunteer lay in the left/right lateral recumbent position for the corresponding leg to be scanned, using the cushion specifically designed to keep the leg at the desired angle. A total of 32 3D MRI volumes were obtained each of 256 × 256 × 192 voxels, where each voxel is of size 0.66 mm × 0.66 mm × 0.69 mm.

3. Methods

For the multi-modal volumetric registration of US with MRI, we propose a two-level initialization i.e., pre-initialisation using structural information followed by segmentation based registration in distance space. This, in turn, initialises US for the final step, LC${}^2$ registration. The flow of the proposed algorithm is shown in Figure 2, and the detailed step-wise procedures are discussed in this section.

3.1. Segmentation

For the multi-modal registration we do in this paper, where we have non-corresponding and unequal number of landmarks, we use segmentations to aid the initialization (Appendix A). The choice of the segmentation algorithm and the number of the anatomical structures segmented generally depend on the specific application and target anatomy. These segmentations can easily be obtained by coarse annotations or fully automatic methods.

3.1.1. Segmentation of Femoral Cartilage from US

We adopt deep-learning based segmentation of femoral cartilage from the US volumes. Deep-learning algorithms provide non-iterative and rapid inference, which is beneficial since the US data are intra-operative and requires near-real-time processing for our target application. We trained a Mask R-CNN instance segmentation network [22] to simultaneously segment the medial and lateral femoral cartilage, also achieving identification along with region delineation [23]. The ground truth labels are marked interactively by the authors and validated by an expert clinician (refer Appendix A.1).

3.1.2. Segmentation of Femoral Cartilage from MRI

Automatic dense labeling of MRI subjects is done by using atlas-based segmentation. An open-source SPL Knee Atlas developed by The Surgical Planning Laboratory at Brigham and Women’s Hospital, Harvard Medical School [24], is elastically deformed onto our MRI volumes. The atlas is an MRI derived atlas generated through semi-automated segmentation algorithms and 3D image reconstruction techniques. Deformable registration of the atlas data onto the subject volume is performed using ANTs (Advanced Normalization Tools). Applying diffeomorphisms, the altas based segmentation computes a reversible transformation composed of a global affine transform and a dense vector field denoting the deformation suffered by the atlas to morph onto the subject MRI image [25]. Figure 3 shows the segmentation output achieved after the deformable registration is applied to the atlas.

3.2. Registration

3.2.1. Data Pre-Processing

Both MRI and US volumes are resampled to obtain an isotropic spacing, i.e., a spacing of 0.33 is chosen and set in all directions for both the modalities. Following that, the US volumes are modified (right-to-left changed to left-to-right or vice versa) in such a way so as to match the anatomical directions of MRI volumes. For example, conventions for left leg are, from the inferior apex of the patella to the superior surface of the tibia, i.e., along the direction of patellar ligament in sagittal plane and from medial femoral cartilage to lateral femoral cartilage in the transverse plane. As an effect of resampling, the size of the MRI image is very large compared to that of the US. It is required for both the modalities to have comparable sizes for close registration and hence an appropriate ROI is chosen from the MRI volume.

3.2.2. Pre-Initialization

During the MIS, US probe is positioned on the patellar tendon and fat pad as shown in Figure 1. We model the US view in the MRI as an oriented bounding box of the anterior knee anatomy. We use the fat pad structure to arrive at a cuboidal region that is close to the view obtained by the US, for the initialization. The infrapatellar fat pad segmentation is obtained by using the method in Section 3.1.2. The principal axes of the fat pad, which are the normals to the faces of the bounding box, are computed from the segmentation, using principal component analysis (PCA) on the points on the fat pad surface.

We now use anatomical information to align the US beam direction to one of the axes. Given the probe placement on the patellar tendon, the fat pad, which is immediately beneath, has almost similar orientation. This enables us to align the inferio-superior direction to one of the fat pad axes directly. The centroid, C, of the fat pad, projected outward gives a point on the anterior fat pad surface, which is our initialization for probe center. Of the three eigenvectors obtained, the vector, $v_{si}$, which goes from superior to inferior direction of the fat pad, gives the rotation matrix for the US volume. Let $p_1$ and $p_2$ be the boundary points of the fat pad that are intercepted by $v_{si}$. The anterior point, $p_1$, gives the required translation for the US volume (refer Appendix B).

Figure 4 shows the fat pad segmentation, the boundary points intercepted by $v_{si}$, and the principal axes of fat pad. O is the centre of the cropped MRI image. The rotation matrix given by $v_{si}$ and the translation vector given by $p_1$ together give us the initial transform of the US.

3.2.3. Initialization by Euclidean Distance Transform

Registering the obtained set of label masks (Section 3.1.1 and Section 3.1.2) extracted from US and MRI volumes can be posed as a pseudo/mono-modal registration problem for which there are many traditional registration approaches available. However, these methods tend to face issues like incorrect error estimation and underperformance due to a lack of efficient optimization. Since labels are encoded with numerical values, one incorrectly registered voxel with label difference of two gives the same error as two incorrectly registered voxels with label difference of one leading to incorrect estimation of registration error. Metrics like Dice are good for finding an overlap of segmentations but might not be sufficient as a metric for registration, since there are equivalences of values at multiple shifts, which would imply multiple local minima in the error surface, which could trap iterative methods. In addition, in this case, a bad initialization might lead to low capture range as the homogeneous label regions such as the background label, do not produce meaningful information for the optimization and therefore result in poor registration [17]. Furthermore, the structures in MRI are complete, whereas, in US, it is dependent on the placement of the probe, so complete volume overlap of the two masks is not expectable from a practical standpoint. Considering these issues, in our method, we transform the masks to a geodesic, using distance transform which makes the support region of the mask larger, thereby increasing capture range. Hence, this enables working with partial segmentations. Each label map is considered as a different binary image, $I_n$, which has two classes foreground $f_g$ and background $b_g$.

Distance transform assigns a number to a voxel which is the distance from that voxel to the nearest non-zero voxel. Euclidean distance transform $\varphi$ is applied on the pre-initialized $I_n$ in Equation (1),

$${\varphi }_{m,n}=\varphi \left(I_{m,n}\right)and{\varphi }_{f,n}=\varphi \left(I_{f,n}\right)$$

(1)

where ${\varphi }_{m,n}$ and ${\varphi }_{f,n}$ are distance maps of moving and fixed modalities respectively, and are derived using Equation (2),

$$\varphi \left(I_n\left(X\right)\right)=\left\{\begin{array}{c}0,I_n\left(X\right)\in f_{g,n}\hfill \\ min(\parallel X-X_0\parallel ),\forall I_n\left(X_0\right)\in f_{g,n}and{I}_n\left(X\right)\in b_{g,n}\hfill \end{array}\right.$$

(2)

$\parallel X-X_0\parallel$ is the Euclidean distance between two points. This transform is applied on all the label classes separately. The initial transform is obtained by minimizing the following Equation (3):

$${min}_{T\in SE\left(3\right)}\sum _{n=1}^N{\int }_{S_f}{|({\varphi }_{m,n}.T)\left(x\right)-{\varphi }_{l,n}\left(x\right)|}^{2}dx$$

(3)

where $T\in SE\left(3\right)$ is the required rigid transformation to initialize ${\varphi }_{m,n}$, $S_f$ is the set of coarse segmentations in fixed volume and N is the total number of segmentations. Since the distance transform of a label map is also dependent on the distance between the segmentation and boundary, the ROI of MRI should be selected accordingly for the minimization problem to work more efficiently. Intuitively, the dimensions and centre for cropping the ROI are selected based on the length and width of femoral cartilage, $fc$. The centre of the segmented femoral cartilage is set as the centre for ROI, and the dimensions are determined by Equation (4):

$$X_{MRI}=\frac{lengthoff{c}_{MRI}}{lengthoff{c}_{US}}·X_{US}$$

(4)

where $X_{MRI}$ denotes the required $x,y,z$ dimensions of ROI, and $X_{US}$ denotes the $x,y,z$ dimensions of US volume and length of femoral cartilage is the length along that particular axis.

3.2.4. Registration

For the registration of the intra-operative US volumes to pre-operative MRI volumes, we use LC${}^2$, a multi-modal similarity metric [11,16]. It is a powerful measure which is invariant to local changes such as brightness and contrast and is proven to work for US-MRI registration [11]. LC${}^2$ correlates both the MRI intensity values and the spatial gradient magnitude with the US. We use this cost function along with optimal patch size to perform the registration. During the registration, it is necessary to look out for and eliminate the unwanted parameter updates in order to avoid unstable results. Therefore, we optimize all the parameters of the obtained rigid transform after every evaluation. It is difficult to compute the analytic derivative of LC${}^2$ which is of least square fitting and hence, as in [11], we use Bound Optimization by Quadratic Approximation (BOBYQA) [26] optimization, which is capable of internally creating its own derivative approximations resulting in fewer evaluations than any other direct search methods.

4. Results and Discussion

To evaluate the performance of the proposed algorithm, an experienced orthopaedic surgeon F.S., manually registered the MRI-US pairs used in this study. The registrations were performed in ImFusion software. The MRI-US volumes were rotated/translated in a common reference system by selecting the rotation/translation values associated with each volume along three orthogonal planes. The software allowed users to visualize the two volumes simultaneously directly while the rotation/translation values were modified by the user.

The MRI volume is kept fixed, while the US volume is rigidly moved to overlay the tibia and the femur surfaces in the two modalities along the three orthogonal planes. The transformation obtained as a result of this process which the surgeon (F.S.) performed manually is considered as the ground truth transform. For every registration, the three translation (mm) and rotation (degree) parameters were obtained from the transform. A similar kind of transform is obtained as a result of registering the two modalities using the proposed algorithm. In-order-to calculate the registration error, the MRI volume is fixed in the position selected in the manual registration performed by the surgeon and the transformation parameters of the US volumes (translation and rotation) are compared as in [27]. The norm of the difference vector of the parameters, computed using Equations (5) and (6), give the translation (t) and rotation (r) error values for each MRI-US pair,

$$t=\sqrt{{(x_{GT}-x_o)}^2+{(y_{GT}-y_o)}^2+{(z_{GT}-z_o)}^2}$$

(5)

$$r=\sqrt{{({\theta }_{GT}-{\theta }_o)}_{yz}^2+{({\theta }_{GT}-{\theta }_o)}_{xz}^2+{({\theta }_{GT}-{\theta }_o)}_{xy}^2}$$

(6)

where subscript $GT$ refers to the ground truth parameter value, and subscript o refers to the outputted parameter value from the developed registration method. x, y, z represent translation parameters in each plane and $\theta$ represents rotation parameter for $yz$, $xz$, $xy$ planes.

The proposed algorithm is tested on 11 MRI-US pairs obtained from six volunteers at four different knee flexion angles. The transform obtained at the end of each of the three stages of the algorithm is compared against the ground truth transform, and the translation/rotation error is computed using Equations (5) and (6).

Table 2. Table showing a decreasing trend in the translation and rotation error (left to right) for each subject computed after each step of the proposed registration algorithm. In-Err: Initial error, Pre-in: Error after pre-initialization, DT: Error after applying Distance Transform, LC${}^2$: error after final LC${}^2$ registration.

Angle	Patient No.	Translation Error (mm)				Rotation Error (${}^{\circ }$)				Final Error(t,r)
		In-Err	Pre-in	DT	LC${}^2$	In-Err	Pre-in	DT	LC${}^2$
0${}^{\circ }$	Subject 1	34.72	9.83	6.04	2.00	73.70	20.13	28.67	5.50	2.00 mm, 5.50${}^{\circ }$
0${}^{\circ }$	Subject 2	35.39	10.59	14.62	3.94	62.14	21.01	14.47	8.78	3.94 mm, 8.78${}^{\circ }$
0${}^{\circ }$	Subject 4	39.10	13.45	2.67	1.72	60.28	8.94	4.93	3.91	1.72 mm, 3.91${}^{\circ }$
0${}^{\circ }$	Subject 5	32.09	7.43	3.88	1.47	65.88	13.29	17.26	4.15	1.47 mm, 4.15${}^{\circ }$
0${}^{\circ }$	Subject 6	31.68	11.04	17.74	2.43	67.83	33.15	13.41	8.67	2.43 mm, 8.67${}^{\circ }$
30${}^{\circ }$	Subject 3	26.14	12.05	1.43	1.41	45.37	32.13	2.41	2.34	1.41 mm, 2.34${}^{\circ }$
30${}^{\circ }$	Subject 6	31.89	30.06	9.19	4.80	25.88	35.29	12.56	9.03	4.80 mm, 9.03${}^{\circ }$
60${}^{\circ }$	Subject 3	27.99	17.75	8.70	5.09	55.21	7.81	15.80	6.43	5.09 mm, 6.43${}^{\circ }$
60${}^{\circ }$	Subject 6	27.71	7.49	7.45	6.00	27.94	18.70	9.70	5.80	6.00 mm, 5.80${}^{\circ }$
90${}^{\circ }$	Subject 3	33.13	22.79	5.36	3.63	25.82	33.29	26.46	15.07	3.63 mm, 15.07${}^{\circ }$
90${}^{\circ }$	Subject 6	34.50	12.15	3.58	6.03	25.09	12.53	17.62	16.04	6.03 mm, 16.04${}^{\circ }$

shows the detailed evaluated error values at each stage from the registration method developed. To calculate the initial error, the transforms of US and MRI volumes are set to identity. It can be clearly seen from Table 2 that there is a decreasing trend in the registration error and that each stage of the registration contributes notably to reduce this error. The pre-initialization step alone, in particular, significantly brings down the initial error by 76.8%, thereby assisting the subsequent steps in the algorithm to achieve global registration more efficiently.

Furthermore, to estimate the importance of including the distance transform step, a cost function analysis is performed. The final cost function, which determines the accuracy of the registration, LC${}^2$, is plotted against any two parameters from the initial transform. We selected the rotation x, rotation z parameters from the output transforms of both pre-in and DT stages, varied their values within the range ±10${}^{\circ }$ and plotted against the LC${}^2$ metric for a given subject. The resulting curves (Figure 5) revealed that the inclusion of distance transform produces a more peaked objective surface (Figure 5b), which would enable better convergence, whereas, without the distance transform (i.e., using LC${}^2$ directly on volumes), the objective surface is undulated with no clear optimum (Figure 5a). In Figure 5b, a global peak is seen with maximum LC${}^2$ value of 0.2, whereas, in Figure 5a, the maximum LC${}^2$ value is observed to be a little lesser with minor local optimum. It can therefore be inferred that incorporating the distance transform step in the algorithm helps achieve a more accurate global registration.

To demonstrate the step-wise decreasing trend in the registration error from Table 2, the registration output at the end of each of the four stages for subject 5 overlaid on their respective ground truths is shown in Figure 6. Figure 6a shows the initial gap between the initial position of the US image and the ground-truth solution provided by the expert and Figure 6b shows the position after the pre-initialization, and it can be seen that this step provides a good initial alignment. The last two steps, Figure 6c,d, help with further fine-tuning and lowering the registration error. It can be spotted from Figure 6d that the patella, the medial femoral cartilage and the patellar shadow of the output red slice overlap almost perfectly with blue ground truth slice. The same registration result is overlaid onto the respective MRI of subject 5 and shown in Figure 7. Furthermore, to assess the efficiency of LC${}^2$, randomized trials were executed on all 11 data sets, repeating the rigid registration for randomly displaced starting positions, i.e., within the range of ±10 mm/${}^{\circ }$ for all six parameters resulting in a total of 1100 registration attempts. Tracing the RMS error and the value of LC${}^2$ metric, we observed that, in many cases, LC${}^2$ is resilient to perturbations in initialization, and hence is an appropriate measure especially for a multi-modal 3D MRI-US registration in a realistic clinical setup. The results for all subjects are shown in Figure 8, including the percentage of converged executions.

To assess the reproducibility of the manual registrations here used as ground truths, a second experienced orthopaedic surgeon (Y.T.${}^5$) performed the manual registrations for a subset of the MRI-US pairs (8 out of 11). The second surgeon utilized the MRI alignment provided by the first expert and then performed the registration between the MRI-US pairs independently from the first surgeon, following the same procedure. Since for each MRI-US pair registered by both surgeons, the MRI volume is fixed in a common position, the difference between the US volume translation/rotation values selected by the two surgeons is utilized to quantify the inter-observer variability as shown in [27]. The inter-observer reproducibility study results reported in

Table 3. Inter-operator variability study.

Angle	Subject	$\mathbf{\Delta }$ Translation	$\mathbf{\Delta }$ Rotation
0${}^{\circ }$	Subject 3	3.00 mm	4.12${}^{\circ }$
0${}^{\circ }$	Subject 5	3.89 mm	6.51${}^{\circ }$
0${}^{\circ }$	Subject 6	4.63 mm	9.24${}^{\circ }$
30${}^{\circ }$	Subject 6	3.64 mm	4.05${}^{\circ }$
60${}^{\circ }$	Subject 3	7.94 mm	9.06${}^{\circ }$
60${}^{\circ }$	Subject 6	1.67 mm	11.99${}^{\circ }$
90${}^{\circ }$	Subject 3	4.77 mm	9.84${}^{\circ }$
90${}^{\circ }$	Subject 6	5.80 mm	7.28${}^{\circ }$

show that the mean difference between the manual registrations performed by the two surgeons is 4.41 mm ± 1.88 mm for translation and 7.77${}^{\circ }$± 2.80${}^{\circ }$ for rotation. The difference in these registrations could be because the rigid match itself might be insufficient for alignment between the knee US and MRI volumes. This is possibly a result of minor leg movements from the reference position during the MRI acquisition. Moreover, during the US acquisition, it is also possible to have small shifts of the patella due to the pressure exerted by the probe that could have possibly affected the overlap of this bone in the two modalities. This effect might have introduced a degree of subjectivity in finding the best compromise to align the bone-like surfaces in the two modalities resulting in a possible mismatch among the operators.

In both our achieved results (Table 2) and inter-operator study (Table 3), it can be observed that the registration error increases as the angle of flexion increases. It can be highlighted from comparing these tables that the registration error of the proposed automated algorithm for almost all tested cases is either smaller or comparable to the inter-observer registration difference (p-value$<0.002$) and thus our registration method can be considered comparable to clinical standards. Figure 9 shows the registration output of subject 6 in frontal, axial and lateral planes for all the four flexion angles. It can clearly be seen from the overlapping images that the registration error (in particular the rotation) keeps increasing as the flexion angle is increased. A possible explanation for this could be, as the flexion angle increases, the lateral and medial femoral cartilage, which are the basis of our algorithm, shift posteriorly and go behind the patella and hence cannot clearly be distinguished in the US volume. In addition, as the knee flexion increases, the area of femoral cartilage in the knee gap decreases, thereby hindering the accuracy of initialization.

5. Conclusions

Automatic registration of pre-operative MRI with intra-operative US enables real-time guidance during the surgery and can enhance the quality and safety of image-guided knee arthroscopy. In this paper, we have introduced a fully automated registration algorithm consisting of a segmentation based pre-initialization, followed by a distance transform step to achieve global initialization. This is followed by registration with an LC${}^2$ similarity metric which uses the total amount of available image information from both the modalities by locally matching US intensities to both MRI intensity and gradient magnitude. Our method is evaluated on 11 MRI-US pairs from six volunteers at different knee flexion angles and resulted in a mean registration error of 3.45 mm and 7.76${}^{\circ }$ in translation and rotation, respectively. In addition, an inter-observer variability study is performed, comparing the registrations performed by two experienced surgeons and resulted in a mean difference between the registrations of 4.41 mm and 7.77${}^{\circ }$. Although the registration error of our algorithm increased with the angle of flexion of the knee, mean registration error of the proposed method is smaller than inter-observer variability (p-value$<0.002$). We have shown that the proposed automated 3D multi-modal registration of US and MRI is globally convergent, robust and efficient. For the US-MRI volumetric registration, the proposed method can form the starting point for developing more complex transforms, including deformable, and can be explored as part of the future work.