Distal Femur Allograft Selection Using a Spectral Shape Descriptor

The automatic selection of bone allografts in a virtual bone databank is an active line of research. This work presents a new approach to solve this problem, based on a recently published image descriptor, called Volumetric Heat Kernel Signature. This descriptor is used to compare the size and shape of three dimensional thresholded computed tomography volumes. This approach is compared to a published method that uses the Iterative Closest Points algorithm to register a segmented search template to different candidates present in the databank. Statistical testing of the agreement between the two methods show that both approaches render similar results in the relevant clinical context. The proposed method avoids incorrect registrations due to local minima and does not require lengthly manual image segmentation and positioning before its use. The new method is conceptually simple and its mathematical basis is sound.


Introduction
In this work we study the problem of automatic selection of allografts of the distal femur in a bone databank composed of thresholded Computed Tomography (CT) images. The aim of this work is to develop an algorithm to automate repetitive measurement tasks executed in long bone tissue banks, as described in a previous paper by Ritacco et al. [1], and to compare it to a previous automatic method published by Sleiman et al. [2].
The measurements obtained in distal fresh frozen femora stored in a donor tissue bank serve as a general size filter for allograft selection. The selection of a replacement bone for a patient under reconstructive surgery is based on different biomechanical compatibility criteria. The first and most basic criterion is the similarity in extent between host and donor bones; in this case the sizes of the bones are evaluated using three linear measurements, each of them between two well-defined anatomical landmarks. These are the transepicondyle distance (A), the anteriorposterior distance in the medial condyle (B), and the anterior-posterior distance in the lateral condyle (C). These landmarks, as shown in Figure 1, can be accurately determined by a trained surgeon from visual inspection of a rough cortical bone threshold based segmentation of a CT scan. The second criterion is the degree of morphological similarity of the bone surfaces. Both criteria and the virtual bone databank workflow are reported by Ritacco et al. [3]. The main aim of this paper is to describe and to evaluate a novel technique that would serve to select an allograft fullfilling both criteria in the context explained above.

Related work
There are a number of different methods for bone databank matching reported in the literature. These methods vary greatly according to the intended application and anatomical area. Paul et al. [4] report a method for hemipelvic allograft selection based on a rigid three dimensional (3D) registration. The registration is performed directly on the CT volumes, using a gradient descent approach to minimize an euclidean metric. The different registrations are validated by visual inspection. Seiler et al. [5] implement a polyaffine registration process to recover six anatomical landmarks on the distal femur that serve as endpoints for the ABC measures. This strategy yields an acceptable estimation of the allograft size, but it does not address the morphology similarity issue. Sleiman et al. [2] describe a surface registration method based on the iterative closest points (ICP) algorithm [6]. The registration error is measured as mean surface distance and as Hausdorff distance. An advantage of this strategy is that both the size and morphology of the bone are taken into account in the error measure. However, this approach has the disadvantage of not being robust enough, since the method used produces registrations where the metrics fall in local minima, but the visual inspection of the result determines an incorrect alignment of the bones. In this work we compare the performance achieved by the mentioned method with the new proposed method.

Contribution
Most of the methods cited above require the existence of an acceptable image segmentation process. This is not always the case since quality segmentation is usually achieved through a long and tedious manual process which is executed by a knowledgeable expert in the domain of anatomical interest. In this work we present as a contribution the application of a spectral based descriptor called Volumetric Heat Kernel Signature (VHKS) to automatically select an adequate distal femur replacement in a bone databank of thresholded CT volumetric scans. VHKS were recently introduced by Raviv et al. [7] and yields useful properties from the viewpoint of a feature descriptor. VHKS are naturally invariant to isometric transforms of a volume. It has been shown in practice that VHKS are also robust to mild non isometric deformations as those commonly found in natural geometric shapes. This is usually the case for distal femora shape variations between subjects. Since VHKS is not scale invariant, it can be used to evaluate both the size and morphology criteria.

Experimental Data and Preprocessing
The experimental data is composed of ten left femora CT volumetric images scanned on a Toshiba Aquilion CT scanner, with a resolution of 0.877 mm and slice increments of 1 mm. The images were thresholded to segment the cortical bone. The segmented images were resampled to a resolution of 5 mm in each dimension, to obtain cubic voxels. The resampled images were cropped to keep only the lower extremity of the femora. Manual segmentations produced by an expert medical surgeon were used to generate surface models, represented as triangular meshes, of each distal femur.

Heat kernel based descriptor
In recent years different spectral methods for shape characterization have been published. Among these methods, the Heat Kernel Signature (HKS) introduced by Sun et al. [8] uses the heat diffusion process on a surface to describe the local and global geometry of each point in the surface with a feature vector. The Volumetric Heat Kernel Signature (VHKS) [7] is derived from the HKS and it is computed on a volume embedded in a three dimensional space. Heat distribution in a domain is governed by the heat diffusion equation, where α is the thermal diffusivity and ∆ X is the Laplace operator. The solution u(x, t) of the heat equation, given adequate initial and boundary conditions, describes the amount of heat on the point x in time t. The solution of the heat equation with an impulse on the surface as the initial condition, u 0 (x) = δ(x − z), is called the heat kernel and can be written as where λ l , Φ l are the eigenvalues and eigenfunctions of the Laplace operator with Neumann boundary conditions [7]. A shape descriptor, the VHKS, is derived by the application of a spatial restriction to the heat kernel If the VHKS of each voxel x of the volume is sampled logarithmically in time (t = γ τ ), then the discrete function is obtained. A vector h x with the first n samples of h x,τ is used as a n-dimensional voxel descriptor, where τ 1 , ..., τ n are different time-scales. The whole volume is described as a normalized sum of the descriptors being #X the number of voxels in X.
The VHKS is invariant to isometric transforms and it is multi-scale, capturing both local features at small t values and global features at large t values. The VHKS is also an informative descriptor, that is, under mild conditions, if two volumes have equal volumetric heat kernel signatures, then these volumes are isometric. However, the VHKS is sensitive to volume scale [8].
The VHKS is computed using a finite differences scheme to evaluate the parabolic partial differential equation that describes the distribution of heat in a body. For each image in the experimental set, the VHKS is calculated for each voxel, starting with an impulse initial condition. Since the heat kernel solution decreases exponentially, each component is normalized against its sum along all the VHKS; also the signature is sampled uniformly at logarithmic time. In this way, the local and global components of the VHKS contribute equally.

Search algorithm
The search algorithm accepts as inputs the VHKS descriptor vector for the template distal femur and the list of descriptor vectors for the bone databank. In this context, the euclidean distance between vectors h X and h Y is used as an index of size and morphology difference between the volumes X and Y . The distal femur volume in the databank whose VHKS vector minimizes this euclidean distance is the best match for the template that is being looked for.

Iterative closest points
The method proposed in this work is compared to the iterative closest points (ICP) registration process. The mean surface distance between two registered surfaces is the metric used to rank the possible matching candidates. Before executing the ICP registration, most of the surfaces were manually positioned to avoid, to a certain extent, the bad alignments due to local minima.   Figure 3: Plots of h X for different volumes X, representing matches found in the ranking generated by the search process for template number 1.
Without this preprocessing step, an elevated number of cases would have been registered in an incorrect way.

Experimental scheme
The agreement between the two methods is measured for the whole experimental set using a leave-one-out scheme: (i) Choose a case to use as template, and leave it out from the experimental set.
(ii) Choose a case to use as candidate, and exclude it from the bone databank.
(iii) Calculate the euclidean distance between VHKS template vector and candidate vector (for ICP, calculate the mean surface distance between the registered surfaces). Store this distance. (iv) If there are still volumes in the experimental set, go back to (ii).
(v) Use the stored distance to rank candidates, from lower to higher magnitudes. (vi) If there are still volumes not used as template, then restore the bone databank and go back to (i).
The process described above produces a total of 100 comparisons and 10 rankings. This scheme is applied to both the VHKS and ICP methods. These rankings are shown in two matrices. Each matrix row is a ranking for certain template; its cells are color coded according to the selection order of the candidate in that column ( Figure 2).

Statistical analysis
Cohen's weighted kappa coefficient is used to measure the inter-rater agreement between the two evaluated methods [9]. The weights are linear and increase indicating the seriousness of the disagreement. For example, weights in the main diagonal of the table of occurrences are 0, weights in the sub and super-diagonal are 1, weights in the diagonals right below and above the sub and super-diagonal are 2, and so on. The statististical analysis and the implementation of the Cohen's weighted kappa coefficient were carried out in the R environment [10].

Results
Cohen's weighted kappa coefficient is κ = 0.45, 95% CI [0.22, 0.68], a moderate agreement according to Landis et al. [11]. If the occurrence table is restricted to the three best matches, then the coefficient is κ = 0.77, 95% CI [0.34, 1.21], a substantial agreement. Color coded matrices showing the different rankings for both evaluated methods are shown in Figure 2. Each row represents a search process. Each column represents a match candidate. The cell color shows the selection order proposed for the template case in its row and the candidate case in its column.
A sample of fours VHKS vectors showing the template case 1, the two best matches and the worst match is shown in Figure 3.

Discussion
The new method proposed in this paper for distal femur automatic allograft selection, produces acceptable results in a short list of best matches. This method is compared to a published method based in the ICP algorithm. The importance of this comparison between ICP and VHKS lies in the fact that the second method is an intrinsic measure of the shape of a volume, whose meaning is not as intuitive as the closest points distances used by the ICP algorithm.
The agreement of the two methods, to a certain point in the application domain, is a step forward accepting the VHKS as a reliable method to select allografts of appropriate size and morphology. Although the general agreement between both methods is classified as "moderate", in terms of the statistical analysis, however the medical relevance of both methods is in the first two or three best matches. These are the selection possibilities the physician is presented with; given this short list he or she will choose the best allograft according to his or her medical expertise. In this scenario the new method may help to obtain the best selection, when restricted to the three best matches, both methods present a "substantial" agreement.
Future developments will refine two aspects of the proposed method. The first one, and most important, is a manual validation of the rankings. In this work the comparison with the ICP based method serves only as a bronze standard. The second aspect is that the low resolution of the grid generated to compute the VHKS using finite differences may show discrepancies in the volume shape that would not be so with finer resolutions. A parallel programming implementation may be able to compute larger problems with a fine-grained grid and with a better performance.
A suitable application of this method may be to select bone allografts in other anatomical regions of the body, like the pelvis. This line of work will be continued in the future.