Improved Surface-Based Registration of CT and Intraoperative 3D Ultrasound of Bones

The intraoperative registration of preoperative CT volumes is a key process of most computer-assisted orthopedic surgery (CAOS) systems. In this work, is reported a new method for automatic registration of long bones, based on the segmentation of the bone cortical in intraoperative 3D ultrasound images. A bone classifier was developed based on features, obtained from the principal component analysis of the Hessian matrix, of every voxel in an intraoperative ultrasound volume. 3D freehand ultrasound was used for the acquisition of the intraoperative ultrasound volumes. Corresponding bone surface segmentations in ultrasound and preoperative CT imaging were used for the intraoperative registration. Validation on a phantom of the tibia produced encouraging results, with a maximum mean segmentation error of 0.34⁡mm (SD=0.26⁡mm) and a registration accuracy error of 0.64⁡mm (SD=0.49⁡mm).


Introduction
e use of image-guided surgery (IGS) systems has expanded significantly for more than 20 years, in the number of procedures performed and the variety of clinical applications, in part due to the continuous advancement of computer power and medical imaging systems (as well as tracking systems, registration methods, software development frameworks, and visualization techniques). e main clinical applications of IGS systems in orthopedics are pedicle screw insertion, hip replacement, knee replacement, and fracture alignment [1].
Scheufler et al. [2] reported the results of image-guided instrumentation of the cervical, upper, and midthoracic pedicles. e study considers the insertion of 248 pedicle screws and reports safe and highly accurate results. Jeswani et al. [3] reported the evaluation of CT image-guided navigation of pedicle screws in small thoracic pedicles. e results showed that image-guided navigation allows for safe, effective, and accurate instrumentation of small ( ≤3 mm) to very small ( ≤2 mm) thoracic pedicles. Deep et al. [4] reviewed the literature on computer-assisted orthopedic surgery (CAOS) for knee and hip arthroplasty; the authors conclude that CAOS systems have gained a pivotal role in lower limb arthroplasty. e use of image-guided surgery (IGS) systems in trauma has been explored for the last eight years approximately, and several methods and systems have been developed with various levels of success [5]. e recent development of minimally invasive arthroscopic procedures for fracture fixation in long bones has substantially improved the accuracy of the procedures with shorter patient recovery times. e success of the surgery depends heavily on the skill and experience of the surgeon [6]; however, the use of CAOS systems in fracture fixation performed with arthroscopic procedures offers an improved visualization of the surgical site with position feedback [7]. It has also been reported that the use of CAOS, in arthroscopic fracture reduction procedures, improves the accuracy of the surgery and achieves a shorter and more effective patient recovery [8]. Buschbaum et al. [9] reported the development of a virtual environment which helps the surgeon with planning optimal reduction paths for femoral fractures. Weil et al. [10] reviewed the literature on the evolution of imageguided iliosacral screw placement, for the reduction of pelvic ring fractures. 3D image-guided surgery is faster and more accurate and uses less X-rays than conventional fluoroscopyguided surgery. e authors stress the challenges of platform interoperability, learning times, and cost reductions. e main stages of most CAOS systems are preoperative image acquisition (usually computed tomography) of the anatomic region of interest, computer-assisted modeling and surgery planning, intraoperative image registration, and finally surgery execution with some computer assistance based on the surgical plan. Intraoperative registration is a critical process of most CAOS systems, given that registration accuracy determines the precision of the surgical visualization, planning, and navigation with respect to the patient on the operating table [11].
Several methods have been developed for intraoperative registration of preoperative CT volumes for CAOS, with fiducial-based registration as the most widely used. Preoperatively, fiducial points are annotated in the medical images or graphic models of the bones, and intraoperatively, the same points are located in the bones of the patient using a navigated tool [12]. Usually only a few anatomical landmarks can be reliably selected; therefore, artificial fiducial markers have been used. e marker is surgically attached directly to the bone of the patient before the preoperative study is acquired. is is an accurate registration approach; however, it increases the complexity of the surgical procedure and the risk of complications due to the invasiveness required to reach each fiducial point.
In order to minimize the invasiveness of the registration process while maintaining high accuracy, image-based methods have been developed. Image-based registration has clear advantages: no need for artificial fiducials, therefore, limiting bone exposure; high accuracy can be obtained through registration of large surface areas. Fluoroscopy-based registration has been extensively used [7,13,14]. However, fluoroscopes are large and expensive and expose the patient and the surgical staff to ionizing radiation. On the other hand, ultrasound is a convenient and safe intraoperative imaging modality; it is portable, low cost, and real time. Its main limitations are low signal to noise ratio due to speckle, user-dependent acquisition and interpretation, and inability to penetrate bones.
Previous studies have shown that the use of intraoperative ultrasound is feasible to achieve errors that satisfy the accuracy requirements of surgery. Herring et al. [15] reported the use of spine phantoms submerged in water to locate fiducial points in ultrasound images. Submillimetric accuracy was reported for the registration of the fiducial points with the spine phantom. Ionescu et al. [16] reported one of the first registration methods based on the segmentation of the bone cortical in ultrasound images. Yan et al. [17] achieved registration errors of less than 2.0 mm in transpedicular screw insertion in porcine cadaver experiments and also less than 1.0 mm of error in phantom experiments.
Most methods for image-based registration of CT and ultrasound can be classified as intensity-based or surfacebased methods. Surface-based registration requires a process of feature extraction in both modalities, and registration accuracy is then affected by feature extraction. Meanwhile, intensity-based registration methods optimize a similarity (objective) function. e following groups, Penney et al. [18], Winter et al. [19], and Gill et al. [20], developed different intensity-based registration methods and reported registration errors under 2.0 mm. e most important difference between both registration approaches (intensity or surfacebased) is its suitability for a specific surgical procedure.
In surface-based registration methods, the manual, or automatic, segmentation of the bone surface in the intraoperative ultrasound images is necessary. Kowal et al. [21] developed an early automatic segmentation method for 2D ultrasound images of bones, based on pixel intensity and position, since bones usually appear bright and are located at the top of the ultrasound images. e method is sensitive to image contrast, having an average segmentation error of 0.42 mm (SD � 0.19 mm).
Beek et al. [22] reported a method for surface-based registration during scaphoid fracture reduction surgery, 3D ultrasound images are segmented semiautomatically with the annotation of 10 seed points in the cortical of the bones, registration is performed using iterative closest points, and a segmentation error of 0.56 mm (SD � 0.08 mm) was reported. Hacihaliloglu et al. [23] developed a bone segmentation method for 3D ultrasound images, based on the Log-Gabor filters and phase coherence, and a segmentation error of 0.28 mm (SD � 0.24 mm) was reported for validation in bone phantoms.
In this work, is reported a new method for the automatic segmentation of the bone cortical in 3D ultrasound imaging. e method is based on a voxel classifier, trained with features obtained from the eigenanalysis of the Hessian matrix corresponding to every voxel in the ultrasound volume. e segmentation of the bone surface in an ultrasound volume is subsequently used for intraoperative registration of the corresponding surface segmented in a CT volume. In the following sections, are reported the acquisition of 3D freehand ultrasound volumes of long bones, the segmentation of the bone surface (cortical), and the registration of a preoperative CT volume. Experimental validation was performed on a phantom of the tibia.

Freehand Acquisition and
Reconstruction of 3D Ultrasound Volumes. 3D freehand ultrasound images were acquired with a conventional 2D scanner (Aloka 1000, Hitachi Aloka Medical America, Inc.) using a 7.5 MHz probe. An optical tracker (Polaris Spectra, Northern Digital, Inc.) was used to navigate the ultrasound probe. Ultrasound B-mode images were continuously acquired using a frame grabber (Epiphan Systems Inc.). For the reconstruction of the 3D ultrasound volumes, the probe was calibrated using the crosswire methodology [24], and 3D volumes were reconstructed using the pixel-based method reported in [25]. Full details of the acquisition and reconstruction of 3D freehand ultrasound images have been published elsewhere [26,27].

Segmentation of the Bone Surface in 3D
Ultrasound. e segmentation method has two main stages: bone surface enhancement, using feature extraction based on differential geometry, and surface segmentation using a Bayes classifier and a region growing algorithm. Bone surfaces in an ultrasound image are identified as bright regions above dark regions caused by the acoustic shadows produced by the bone. High-intensity echoes produced by the bone surface can be observed in ultrasound images as bright lines with a width of 2-4 mm. Jain and Taylor [28] have shown that a good estimate of the bone cortical in an ultrasound image lies at the center line of the bright echo line. e center line of the ultrasound echoes corresponds to the maximum principal curvature in the, orthogonal, gradient direction.
is allows for the detection of the bone cortical using a method based on ridge detection as described below.

Bone Surface Enhancement in 3D Ultrasound.
e second derivatives, in each direction around a voxel in an ultrasound volume, provide information about the type of geometric shape to which the voxel belongs. e geometric shapes could be blobs, tubes, or surfaces [29]. Shape information was obtained from the principal component analysis of the Hessian matrix calculated on each voxel I(x, y, z) in the volume. e Hessian matrix (2) was constructed with the second partial directional derivatives, which are approximated with the convolution of every voxel in the volume I(x, y, z) with a Gaussian kernel (1) in each direction, where σ is the scale factor of the kernel. e size of the kernel is determined as three times σ. e Hessian matrix of a voxel is a symmetric matrix that contains the partial derivatives in all the possible directions as shown in the following equation: where each z ij I represents the second order partial derivative of the voxel in the directions i and j. e principal component analysis of the Hessian matrix results in three eigenvalues λ 1 , λ 2 , and λ 3 and the corresponding eigenvectors v 1 , v 2 , and v 3 , for each voxel in the ultrasound volume. If λ 1 ≥ λ 2 ≥ λ 3 , there is a set of conditions that determines the membership of a voxel to a certain type of geometric shape which can be tubular, spherical, or a surface. e conditions are shown in Table 1.
Based on the conditions shown in Table 1, voxels that belong to a surface were chosen and assigned a new value given by (3). is highlighted all the voxels that had a high probability of belonging to a surface [29], Equations (3) and (4) represent the condition that one voxel belongs to a surface. e condition λ 3 ≪ λ 2 ≃ λ 1 ≃ 0 was expressed as λ 3 ≪ λ 2 ≃ 0 and λ 3 ≪ λ 1 ≃ 0; those inequalities are described by (4) where the value of ω decreases with the deviation from the condition λ 1 ≃ 0. As recommended by Sato et al. [29], we choose the values α � 0.5 and c � 1.

Bone Surface Extraction in 3D Ultrasound.
A bone segmentation algorithm based on a Bayes classifier was developed. e classifier was trained with features of the voxel intensity and the first, second, and third moments of the enhanced voxel values previously calculated. ree classes were considered for training: bone, soft tissue, and acoustic shadow produced by the bone. For each voxel, a feature vector was constructed as follows: where I(x, y, z) and E(x, y, z) are the original image and the enhanced image values at a specific position; μ, σ 2 , and c are the mean, variance, and skewness in a (9 × 9 × 9) window centered in E(x, y, z). For each vector x → (x, y, z), the probability of membership to each class, bone, tissue, and shadow, was estimated using the Bayes posterior probability defined in the following equation: where P(C k ) is the a priori probability for each class obtained from the training sample proportions and P( x → ) is the a priori probability of vector , we assumed a normal distribution for each class (bone, tissue, and shadow) as defined in the following equation [30]: where is the covariance matrix and μ → is the mean of a set of vectors x → that belong to the class. Substituting (6) in (5) and discarding P( x → ), we can calculate a discriminant function Y k : Equation (7) allowed us to classify each voxel in an ultrasound volume as follows: for each voxel, we calculate the corresponding vector x → and evaluate (7) for each of the classes: tissue, bone, and shadow. A voxel was labeled as bone (i.e., labeled as " ; otherwise, it was labeled as background (i.e., labeled as "0"). e mean and covariance of each class (bone, tissue, and shadow) were estimated on a small training volume (150 × 150 × 150) which should include voxels of the three classes. To produce one complete surface with a minimum of holes, a region growing process was also applied [31]. Seed points were selected from the 0.01 percent of the previously classified bone voxels with the highest probability. Starting from these seed points, the bone surface was segmented considering all the neighboring voxels previously classified as bone.

Intraoperative Registration of CT Volumes.
e bone surface (cortical) was segmented manually in the preoperative CT volume, and automatically in the intraoperative ultrasound volume, using the method reported in the previous section. From each segmentation, the corresponding mesh was generated using Marching Cubes [32], and each mesh has around 2.5 × 10 5 vertices. Registration was performed in two stages. First, a few points (three or four) are selected manually with a coarse accuracy, in both meshes, to have an initial rough alignment required by the iterative closest point (ICP) method. en a subset of 2000 points was randomly acquired from both meshes. Mesh registration was performed using iterative closest points [33]. e software was developed using 3D Slicer [34] with specific modules developed for this research using ITK (http:// www.itk.org) and VTK (http://www.vtk.org).

Results and Discussion
e accuracy of the bone segmentation and registration was evaluated on a realistic phantom of the tibia, which was constructed with a synthetic tibial bone (ERP #1117-42, Sawbones Inc., Vashon, WA, USA) immersed in a hydrogel made of polyvinyl alcohol (PVA) diluted in 95% of water. It has been shown that PVA resembles the appearance and the mechanical properties of soft tissue on ultrasound images [35]. A passive tracking tool was firmly attached to the distal end of the tibia phantom. Figure 1(a) shows a photo of the PVA phantom of the tibia. e dimensions and shape of the tibia phantom were accurately measured using microCT.

CT Image Acquisition.
MicroCT (Nikon Metrology XTH 225) with a 2048 × 2048 pixel matrix was used to scan the phantom of the tibia including the tracking tool; imaging settings were 220 kV, 61 µA, 3142 projections, and one image per projection, and no physical filter was used. e final volume dimensions were 1042 × 1250 × 3201 voxels in (x, y, z) axes, respectively, with an isometric voxel size of 0.115 mm. e acquired CT volume had a very high resolution that helps improving the accuracy measurements due to a small resolution error (0.057 mm) in the CT data. Figure 1(c) shows three orthogonal views of the acquired CT volume of the phantom.

3D
Freehand Ultrasound Image Acquisition. 3D freehand ultrasound images were acquired from the tibia phantom, and the volume was reconstructed as described in Section 2.1. e origin of our coordinate frame was located at the center of one of the reflecting spheres of the tracking tool attached to the tibial bone, shown in Figure 1(a), and the voxel size for the reconstruction of the ultrasound volume was the same of the CT: 0.115 mm.
e phantom was scanned with the ultrasound probe using the same waterbased gel used for clinical ultrasound scans. e phantom was scanned at different sections always including the tibial plateau. In Figure 1(b), are shown three orthogonal views of the acquired ultrasound volume of the phantom.

3D Ultrasound Segmentation.
A small ultrasound volume (150 × 150 × 150 voxels) was acquired from the test phantom, for training of the segmentation method. Bone, tissue, and shadow regions were manually annotated. e corresponding voxel feature values were stored as the training sample for the Bayes classifier described in Section 2.2. is training set was used for all the experiments reported below, and the training process was performed only once for all experiments. e accuracy of the bone segmentation of the intraoperative ultrasound volumes was measured using as a reference the manual segmentation of the bone in the high-resolution CT data. is segmentation was approved by an expert orthopedic surgeon. As previously mentioned, all the intraoperative ultrasound images were acquired with the origin of the navigation reference frame located at the tracking tool attached to the plastic bone. Since the position of the tracking spheres in the corresponding CT volume was determined, the CT and the ultrasound images were accurately registered. All the remaining errors between the bone surface in CT and ultrasound were mainly due to segmentation errors of the bone cortical in the ultrasound volume (as illustrated in Figure 2).
Twelve acquisitions of 3D freehand ultrasound volumes were performed at different sections of the phantom of the tibia, which was located each time at different positions and orientations on the experiment table. e bone cortical on each volume was automatically segmented, and the corresponding mesh was constructed. Table 2 summarizes the mean distances of each node of the ultrasound mesh to the nearest node of the reference CT mesh. As can be observed in Table 2, the maximum mean segmentation error was 0.334 mm (SD � 0.257 mm) for experiment four, and Figure 3 shows graphically the segmentation errors. e overall mean segmentation error for the twelve experiments performed was 0.21 mm (SD � 0.061 mm). e average processing time per experiment including volume reconstruction and volume segmentation was about five minutes. Holes in the segmentation, due to shadows in the acquisition of bone ultrasound, were not taken into account.

Intraoperative Registration.
Twelve registration experiments were performed; on each one of them, the mesh obtained from the manual segmentation of the bone surface in the CT images of the phantom was registered with the tibia phantom located at different positions on the experiment table; using the method described in Section 2.3, the errors of the registration process were measured as follows.
For each position of the phantom on the table, a 3D freehand ultrasound volume of a different section of the tibia was acquired. e bone surface in the ultrasound volume was automatically segmented using the method described in Section 2.2.2; from the resulting segmentation, a mesh called M US-intraop was constructed. e mesh constructed from the  (Figure 4(a)). In order to have a reference to measure the registration error, M CT-preop was also registered with the phantom on each position using the transformation given by the optical tracker, and this accurately registered CT mesh was called M CT-ref (Figure 4(b)). e accuracy of each registration experiment was then measured between the registered CT mesh (M CT-intraop ) and the reference CT mesh of the tibia phantom (M CT-ref ), as described below.
For each experiment, the same preoperative CT mesh M CT-preop was registered against the tibia phantom by two independent methods: M CT-intraop registered with the method reported in Section 2.3 and M CT-ref exactly registered using the tracking tool. Target registration errors (TREs) of each experiment were measured as the mean distance between all corresponding vertices in both registrations. Figure 5 shows the results of one experiment before ( Figure 5(a)) and after ( Figure 5(b)) registration, the reference mesh is shown in yellow, the ultrasound generated mesh is shown in blue, and both meshes are shown overlapping on one slice of the ultrasound volume. Figure 6 shows the errors for six experiments, showing the ultrasound   Table 3 and are shown graphically in Figure 7. e average time taken by the registration process was approximately two minutes.

Discussion.
A fully automatic 3D method was developed for the segmentation of the cortical bone in 3D ultrasound images, acquired with the 3D ultrasound freehand technique, which enables the acquisition of large volumes. e eigenvalues, corresponding to the principal component vectors, of the 3D Hessian matrix of each voxel in an ultrasound volume, were used to enhance the bone surface following the work of Sato et al. [29]. A Bayes classifier was trained with five features: the original and enhanced voxel intensity values, as well as the first, second, and third moments of the enhanced voxel values. Validation of the segmentation method on a realistic phantom of the tibia immersed in PVA, allowed for accurate estimates of the segmentation errors, since PVA has mechanical properties similar to those of tissue in ultrasound imaging. e acquisition of ultrasound images was performed with the same water-based gel used clinically.
Segmentation errors were measured as the distances between all the points of the segmented ultrasound mesh and the nearest points in the reference CT mesh (Section 3.3). Table 2 summarizes the segmentation errors obtained for twelve different ultrasound volumes of the same phantom (approximate volume size of 850 × 850 × 850 voxels that corresponds to 10 × 10 × 10 cm). A maximum mean error of 0.334 mm (SD � 0.257 mm) for experiment four is shown. Figure 3 shows the distribution of all the mean error values, and the overall mean for the twelve segmentation experiments performed was 0.21 mm (SD � 0.061 mm). ese results show improvement over previously reported methods for the segmentation of bones in 3D ultrasound. Kowal et al. [21] reported the automatic segmentation of bone contours in 2D ultrasound images, with a mean segmentation error of the bone contour lines of 0.42 mm (SD � 0.19 mm). Beek et al. [22] reported an IGS system for the fixation of scaphoid fractures based on preoperative CT and intraoperative ultrasound images. A semiautomatic method was developed for the segmentation of bone contours in 2D ultrasound. Validation was performed on plastic phantoms of the scaphoid bone immersed in water. e authors reported a mean segmentation error of 0.56 mm (SD � 0.08 mm). Hacihaliloglu et al. [23] reported an automatic method for bone segmentation based on Log-Gabor filters and phase coherence. A mean segmentation error of 0.28 mm (SD � 0.24 mm) was reported for validation on a realistic bovine bone phantom immersed in solid gel. e bone surfaces, obtained from the automatic segmentation of the ultrasound volumes, were used to register a highresolution CT model of the tibia using iterative closest points. Twelve intraoperative registration experiments were performed using the same tibia phantom located at different positions on the experiment table. e target registration errors (TREs) of each experiment were calculated as the distances between all corresponding vertices in the registered CT mesh (M CT-intraop ) and a reference CT mesh (M ct-ref ). Table 3 shows a maximum mean TRE of 1.541 mm (SD � 0.229 mm) for experiment two. e overall mean TRE for the twelve experiments was 0.64 mm which is smaller than other TREs previously reported: in the work of Beek et al. [22], was reported a mean TRE of 3.32 mm (SD � 1.41 mm) for three alumina beads of 3.0 mm of diameter, which includes the uncertainty to locate the exact center of the alumina beads; Penney et al. [18] reported a RMS  Figure 6 shows the target registration errors (TREs) illustrated in a false color scale, for each point on the surface of the CT. As expected, smaller registration errors are obtained when large areas of the tibia phantom are scanned with the ultrasound probe. e lower row of Figure 6 shows three cases with mean target registration errors smaller than 0.24 mm. e shape of the diaphysis (i.e., the middle section) in long bones is very similar along the bone. is makes it difficult to achieve an accurate registration if the scanned area contains only the diaphysis of the bone. In order to obtain better accuracy in the registration process, it is advisable to include part of the epiphysis (i.e., ends of the long bone) into the scan since the epiphysis contains features that can be captured with ultrasound images. Figure 6 shows that when the scanned areas do not contain part of the epiphysis, the registration error is larger. e mean registration time including ultrasound volume reconstruction, surface segmentation, and ICP registration was approximately eight minutes. All computations were made on a Mac Pro with a 2.8 GHz Quad Core Intel Xeon processor and 16 GB of RAM.

Conclusions
Fast and accurate intraoperative registration is a critical stage of most computer-assisted orthopedic surgery (CAOS) systems. e use of intraoperative ultrasound for image-guided registration in CAOS has several advantages: low cost, no exposure  to ionizing radiation, and compact size. However, the accurate registration of CT and ultrasound of bones is still a research challenge due to the low signal to noise ratio of ultrasound imaging and its incapability to penetrate bones. A new method for the intraoperative registration of preoperative CT volumes of long bones was reported in this work. e method is based on the automatic 3D segmentation of the bone cortical in 3D freehand ultrasound imaging which, in turn, enables the acquisition of large bone sections. e method is able to segment the bone surface in large ultrasound volumes with an overall mean segmentation error of 0.21 mm (SD � 0.061 mm). e approximate segmentation time per volume was 1 min (for volumes of 850 × 850 × 850 voxels). e corresponding preoperative CT was manually segmented, and both meshes (ultrasound and CT) were accurately registered using iterative closest points. Accurate registration was achieved, with minimum user interaction. e maximum TRE was under 1.55 mm, and the mean maximum TRE was 0.64 mm.
It has been observed that long bones have very similar surface shapes in completely different locations. A conventional nonnavigated 3D ultrasound probe can only scan a small part of the bone. is makes 3D freehand ultrasound as an optimal choice for intraoperative registration of long bones since large parts of the bone can be scanned and automatically segmented with the method reported. After the acquisition of the 3D ultrasound images, the total processing time was approximately 8 min, with 5 min. for volume reconstruction, 1 min. for bone segmentation, and 2 min. for CT registration. Ultrasound volume reconstruction is suitable for parallel implementations which can reduce significantly the total reconstruction times. Atesok et al. [5] reported that, on average, 14 min. of extra surgical time are added for 2D fluoroscopic navigation in fracture reduction surgery.
is extra time is acceptable for the surgeons given the localization advantages of navigated instruments during minimally invasive procedures. e image segmentation and registration methods reported here are suitable for minimally invasive surgery of long bones such as fracture reduction of the femur and tibia.
Other arthroscopic surgical procedures such as total knee replacement will be subject to the possibilities to scan the surgical site with ultrasound imaging.