PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2017; 12(11): e0187874.
Published online 2017 November 27. doi:  10.1371/journal.pone.0187874
PMCID: PMC5703488

A novel registration-based methodology for prediction of trabecular bone fabric from clinical QCT: A comprehensive analysis

Vimal Chandran, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing,* Mauricio Reyes, Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing, and Philippe Zysset, Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing
Deepak Vashishth, Editor

Abstract

Osteoporosis leads to hip fractures in aging populations and is diagnosed by modern medical imaging techniques such as quantitative computed tomography (QCT). Hip fracture sites involve trabecular bone, whose strength is determined by volume fraction and orientation, known as fabric. However, bone fabric cannot be reliably assessed in clinical QCT images of proximal femur. Accordingly, we propose a novel registration-based estimation of bone fabric designed to preserve tensor properties of bone fabric and to map bone fabric by a global and local decomposition of the gradient of a non-rigid image registration transformation. Furthermore, no comprehensive analysis on the critical components of this methodology has been previously conducted. Hence, the aim of this work was to identify the best registration-based strategy to assign bone fabric to the QCT image of a patient’s proximal femur. The normalized correlation coefficient and curvature-based regularization were used for image-based registration and the Frobenius norm of the stretch tensor of the local gradient was selected to quantify the distance among the proximal femora in the population. Based on this distance, closest, farthest and mean femora with a distinction of sex were chosen as alternative atlases to evaluate their influence on bone fabric prediction. Second, we analyzed different tensor mapping schemes for bone fabric prediction: identity, rotation-only, rotation and stretch tensor. Third, we investigated the use of a population average fabric atlas. A leave one out (LOO) evaluation study was performed with a dual QCT and HR-pQCT database of 36 pairs of human femora. The quality of the fabric prediction was assessed with three metrics, the tensor norm (TN) error, the degree of anisotropy (DA) error and the angular deviation of the principal tensor direction (PTD). The closest femur atlas (CTP) with a full rotation (CR) for fabric mapping delivered the best results with a TN error of 7.3 ± 0.9%, a DA error of 6.6 ± 1.3% and a PTD error of 25 ± 2°. The closest to the population mean femur atlas (MTP) using the same mapping scheme yielded only slightly higher errors than CTP for substantially less computing efforts. The population average fabric atlas yielded substantially higher errors than the MTP with the CR mapping scheme. Accounting for sex did not bring any significant improvements. The identified fabric mapping methodology will be exploited in patient-specific QCT-based finite element analysis of the proximal femur to improve the prediction of hip fracture risk.

Introduction

Osteoporotic hip fractures represent a major clinical and public health problem in aging populations. Identifying individuals at higher fracture risk would enable targeted osteoporosis management and improve fracture prevention. Areal bone mineral density (aBMD) measured by dual-energy x-ray absorptiometry (DXA) is routinely used as a surrogate of bone strength for osteoporosis diagnosis and fracture risk assessment. Modern techniques such as finite element (FE) analysis allow for a more accurate estimation of bone strength using the local distribution of BMD provided by QCT, but do not account for the anisotropy of trabecular bone architecture called fabric. Recent validation studies have demonstrated that the inclusion of bone fabric (anisotropy) in FEA models is important and delivers an improved prediction of bone strength [15]. However, measuring bone fabric requires high resolution peripheral QCT (HRpQCT) images and presently, this resolution is not available clinically for the proximal femur.

Consequently, computational approaches to accurately predict bone fabric directly from clinical QCT images are receiving increasing interest. In this regard, machine learning approaches have been recently used to predict bone fabric, where the statistical relationship between clinical QCT imaging information, and its corresponding high-resolution peripheral QCT (HRpQCT) was modelled, and then used to perform inference of bone fabric on (unseen) clinical QCT images. Particularly, discriminative models that infer bone fabric from computed features (i.e predictor variables) have been proposed. In [6], nodal displacements of a template mesh registered to a patient-specific mesh are used as features for a non-linear kernel partial least square (PLS) regression approach. In [7], morphology- and texture-based features are used as features as part of a decision forest regression approach. While these statistical approaches have showed promising results, they involve manual annotations of landmarks for initial alignment, and their accuracy depends on the selected training data.

Another family of approaches for bone fabric estimation is based on image or mesh registration. In [3], the authors rely on a database of HRpQCT-based derived FE models of femurs including bone density and fabric information. From the database, the most similar femur to the target femur is selected by means of mesh-morphing and a bone-mineral-density similarity metric computed across the database. Finally, the pre-computed fabric information of the selected femur is mapped to the patient’s femur by rigidly correcting the local orientation of the fabric information. It is noted, however, that this study did not perform a direct study on clinical CT images. Recently, in [2, 8], the authors investigated intensity based registration methods to derive fabric information. In their approach, rather than pre-computing fabric information and then mapping the closest femur in the database to the patient’s femur, bone fabric is inferred by registering a single QCT image to the patient’s image, and then it calculates fabric on the corresponding non-rigidly transformed HRpQCT image.

In these previous studies two issues are identified with respect to the chosen registration approach and the degrees of freedom of the transformation model used to map fabric information to the patient’s image. The study in [3] used a surface-based mesh morphing approach (using eight sparsely located landmarks), and a rotation-based local correction to map fabric information to the target image. It is first remarked that surface-based registration approaches have been reported to be less accurate than intensity-based registration approaches for establishing anatomical point correspondences [9]. Secondly, the study in [3] uses a local correction based on a rotation matrix, which is not proved to provide the best result in terms of fabric matching to a patient’s image. Similarly, the approaches in [2, 8] employ an image-intensity registration approach and the complete non-rigid transformation (i.e. no decomposition or local correction of the transformation) to derive fabric information. In this regard, as demonstrated in the present study, the degrees of freedom of the transformation model used to map fabric information to the patient’s image plays an important role on the accuracy of these methods.

Consequently, and differently from previous approaches, we propose a novel registration-based estimation of bone fabric directly from clinical QCT image. It is designed to preserve tensor properties of bone fabric and to map bone fabric by a global and local decomposition of the gradient of a non-rigid image registration transformation. Another issue investigated in this study comes from the fact that the role of utilising a database of femoral atlases, from which fabric information is mapped to a patient’s image, is not known and inconclusive from the state of the art. The conclusions presented in [2, 8] contradict with those of [3], on the fact that a single femur atlas might suffice to estimate femur fabric from a QCT patient image. These contradictory results might be amplified by the fact that these studies were tested on a very limited set of ten cases. In this study we therefore present a thorough leave-one-out analysis on the importance of atlas selection for bone fabric estimation on a dataset comprising 36 pairs of QCT and HRpQCT human femora. Using a deformation-based distance metric, we evaluate six different atlases that span different degrees of similarity to the target image, and are population-wide or sex-specific. In addition, beyond bone shape and image atlases, we evaluate and report on the ability of a single population-, and sex-specific atlases of bone fabric used within the proposed registration-based fabric estimation approach.

Methods

In this section the proposed image registration based fabric prediction is presented, followed by the methodology and metrics proposed to select a femur atlas from a given population. The section continues then with the proposed methodology to decompose the image transformation and apply it to the precomputed fabric of the chosen atlas. The section finishes with the scheme and metrics used to evaluate the quality of the proposed fabric prediction approach.

The complete overview of the registration approach in bone fabric prediction is presented in Fig 1 and is described in detail below.

Fig 1
Overview of the image registration approach for predicting the bone fabric is presented.

Image registration

Image registration is the process of aligning two images into a common coordinate system. Given a pair of images, a fixed image IF(x) and a moving image IM(x) are defined on their own spatial domain: ΩF ⊂ ℝ3 and ΩM ⊂ ℝ3, and here x = {x1, x2, x3} denotes the voxel location. Image registration is the task of finding a coordinate transform T:ℝ3 → ℝ3 that spatially aligns the two images such that a given similarity metric between IF(x) and IM(T(x)) is optimized [10]. Image registration can be formulated as an optimization problem:

T^=argminC(T;IF,IM)=argmin(Csimilarity(T;IF,IM)+γCsmooth(T)).
(1)

The cost function C defines the quality of alignment, which is separated into a similarity measure Csimilarity and a regularization term Csmooth. In this work, normalized correlation coefficient [10] is used as the similarity measure because of its ability to handle mono-modal image registration. Curvature regularization [11] is used as regularization term to cope with the ill-posedness of the non-rigid image registration. It acts on the deformation field computed on the B-Spline grid nodes. The parameter γ weighs regularity against similarity.

Csimilarity(T;IF,IM)=xΩ(IF(x)IF¯)(IM(T(x))IM¯)(IF(x)IF¯)2(IM(T(x))IM¯)2,
(2)
Csmooth(T) =  ∥ ΔxT2
(3)

In the present study the image registration process is performed in two stages. First, an affine registration is performed to get a coarse global alignment of the entire anatomy. Second, a cubic B-Spline registration is used to yield a fine local alignment based on a grid of J control points. The transformations are combined by composition, as follows

T(x) = TB(TA(x))), 
(4)

where TA is the affine transform and TB is the B-spline transform.

Parameter tuning of the registration was performed heuristically and based on the quality of the registration. To this end, we computed the Dice coefficient between image masks, which are obtained via semi-manual segmentation of the HRpQCT images for which a simple image thresholding is feasible. The Dice coefficient is then calculated on image masks transformed (i.e. Eq (4)) and resampled to the QCT image space by nearest-neighbor interpolation. The accuracy of the image registration in terms of Dice coefficient [12, 13] was in average of 94±3%, and hence considered satisfactory for the rest of the analyses. Furthermore, changing the order of operation (i.e. HRpQCT masks were first resampled and then transformed for Dice coefficient calculation) did not significantly affect the accuracy of the transformation (p>0.05).

Selecting a femur atlas

A femur atlas is a QCT image chosen from the population. Since there are various possible candidate femur atlases available in the population, we propose a strategy for choosing it. In principle, a good atlas is such having minimal image deformation needed to warp the atlas image to each fixed image in the population. Inspired from the Frechet mean and related works proposed in computational anatomy [1416], a distance metric DM is proposed herein to measure the extent of deformation. In the proposed image registration process (Fig 1), IF corresponds to the patient’s femur QCT image, while IM corresponds to the femur atlas QCT image.

The distance metric DM is calculated using the stretch tensor jGV, which is computed on the grid of control points {j = 1,…J} spanned over the entire registered image. The computation of the stretch jGV involves the combined transformation of affine and B-spline transforms. We use jGV [equivalent] V for simplicity in the rest of the paper.

The deformation gradient F is computed, which is the gradient of the transformation or the Jacobian matrix of the mapping

F = ∇xT.
(5)

Performing VR decomposition

FVR
(6)

V = (FFT)1/2
(7)

The distance metric DM is defined as

DM=j=1J||log(jV)||=j=1J||log(jVBjVA)||,
(8)

where VA and VB denotes the principal stretch of the affine and B-spline transforms, respectively.

The distance metric DM is then used to select different atlases featuring different degrees of similarity to the target fixed image. Six different and representative femur atlases were chosen to evaluate the importance of selecting an appropriate femur atlas (see Fig 2). For concision, they are henceforth referred to as:

Fig 2
Femur atlas selection strategy.
  1. Closest to the patient femur in the population (CTP,CSP): This atlas image corresponds to the femur image yielding the minimum distance metric (hence referred as closest to the patient’s femur). If N represents the total number of femurs in the population, then the Closest to the patient’s femur in the Total Population, termed here CTP, is
    CTP=min(q=1NDM(patientIF,qIM)).
    (9)
    Similarly, if Nsex represents the total number of femurs in the sex-specific population, then the Closest to the patient’s femur in the Sex-specific Population, termed here CSP, is
    CSP=min(q=1NsexDM(patientIF,qIM)).
    (10)
  2. Farthest to the patient femur in the population (FTP,FSP): The femur image yielding the maximum distance metric is considered to be the farthest to the patient’s femur. If N represents total number of femurs in the population, then the Farthest to the patient’s femur in the Total Population, termed here FTP, is
    FTP=max(q=1NDM(patientIF,qIM)).
    (11)
    Similarly, if Nsex represents the total number of femurs in the sex-specific population, then the Farthest to the patient’s femur in the Sex-specific Population, termed here FSP, is
    FSP=max(q=1NsexDM(patientIF,qIM)).
    (12)
    We note that inclusion of this femur as potential atlas is meant to provide a worst-case scenario, where the atlas and the patient’s femur are considerably geometrically different.
  3. Mean femur of the population (MTP,MSP): Generally, the mean femur of the population is a synthetic image produced through arithmetic computation [17]. However, such synthetic images are prone to present blurred intensity patterns of the femur fabric, stemming from the averaging process. Hence, we chose as mean femur atlas, the real femur image yielding the minimum accumulated distance metric across the population. If N represents the total number of femurs in the population, then the Mean femur of the Total Population, termed here MTP, is
    MTP=min(p=1Nq=1pqNDM(pIF,qIM)).
    (13)
    Similarly, if Nsex represents the total number of femurs in the Sex-specific population, then Mean femur of the Sex-specific Population, termed here MSP, is
    MSP=min(p=1Nsexq=1pqNsexDM(pIF,qIM)).
    (14)

Bone fabric extraction

In this section we briefly describe the step of extracting and modeling fabric information. Bone fabric describes the preferential alignment and structural anisotropy of bone trabecular micro-architecture. It is computed using the MIL method [18], which measures the average distances of bone-marrow interfaces in multiple orientations on a segmented image. In summary, a Laplace Hamming filter is first applied to sharpen the HRpQCT image, which is then normalized, and segmented based on image thresholding [19]. On the segmented image, a cubic volume of interest (VOI) with a side length of 5.3mm is extracted, at each corresponding control point IM(T−1(jx)), from the trabecular region, and fabric tensor jM is computed using the MIL method.

The resulting spatial distribution can be described with a second-order fabric tensor M ∈ ℝ3X3 with eigenvalues mi and normalized eigenvectors mi.

M=i=13mi(mimi)
(15)

where m1m2m3. The fabric tensor M is normalized by dividing it by its trace and multiplying it by a factor of 3 such that

tr(M) = 3.
(16)

The shape of the fabric tensor can be visualized as an ellipsoid with magnitude of eigenvalues providing the indication of the extent to which the structure is preferentially aligned. An elongated ellipsoid represents an anisotropic structure (high degree of anisotropy) whereas a sphere represents an isotropic structure (absence of anisotropy).

Fabric tensor mapping

Computing the fabric tensor directly on the atlas image, which is transformed to the patient’s image via registration (e.g. as in [2, 8]) might result in loss of information as the registration process, involving local image deformations, tends to alter the bone fabric pattern. Contrarily, rather than computing fabric tensors on a transformed atlas image, we propose to map fabric information from the atlas to the patient image by transforming its tensorial representation instead. This is inspired by similar strategies followed in neuroimaging, for DTI image registration [20], where structural MRI is used for an initial registration and then diffusion tensor information is mapped based on the resulting transformation. This is mainly performed to reduce shape variance and to maintain direction consistency.

Consequently, fabric tensor mapping is modeled as the coordinate transform T:ℝ3X3 → ℝ3X3 involved in transforming the fabric tensor M from the space of the femur atlas image to the space of the patient’s femur image. The image registration process involves global and local deformations, which can be decomposed into stretch and rotation components. Understanding the impact of different components of deformations on tensor mapping becomes essential. In this regard, we have chosen five different tensor mapping schemes reflecting different degrees of freedom of the transformation used for fabric tensor mapping. For concision, they are henceforth referred to as:

  1. No Rotation (NR): Fabric tensor mapping involves only translation, which is a direct mapping from the femur atlas to the patient’s image. After image registration, point correspondences are established between patient femur and femur atlas. If M represents the computed fabric tensor from the femur atlas HRpQCT image, then tensor mapping based on No Rotation, termed here NR, is
    MIMITM.
    (17)
    We note that inclusion of this mapping method is meant to show the impact of tensor mapping and its advantages.
  2. Affine Rotation (AR): Fabric tensor mapping involves affine rotation, which is a global transformation. After image registration between the patient image and the atlas image, the affine rotation matrix RA is derived from the deformation gradient F. If M represents the computed fabric tensor from the femur atlas HRpQCT image, then tensor mapping based on Affine Rotation, termed here AR, is
    M=RAMRAT.
    (18)
    We note that tensor mapping by AR will not alter the eigen-values mi but only eigen-vector mi of M.
  3. Affine Deformation (AD): Fabric tensor mapping involves affine deformation, which is a combination of an affine rotation matrix and an affine stretch tensor, and it is a global transformation. After image registration between the patient image and the atlas image, the affine rotation matrix RA and affine stretch tensor VA is derived from the deformation gradient F. Then, the affine deformation gradient FA = VARA is computed. If M represents the computed tensor from the femur atlas HRpQCT image, then tensor mapping based on Affine Deformation, termed here AD, is
    M=FAMFA1.
    (19)
  4. Complete Rotation (CR): Fabric tensor mapping involves complete rotation, which is a combination of an affine rotation matrix and a B-spline rotation matrix. This is a local transformation. After image registration between the patient image and the atlas image, the affine rotation matrix RA and B-spline rotation matrix jRB is derived from the deformation gradient jF. Then, the complete rotation matrix jR = jRB [composite function (small circle)] RA is computed. If jM represents the computed tensor from the femur atlas HRpQCT image, then tensor mapping based on Complete Rotation, termed here CR, is
    jMjRjMjRT.
    (20)
    We note that tensor mapping by CR will not alter the eigen-values mi but only eigen-vector mi of M.
  5. Complete Deformation (CD): Fabric tensor mapping involves complete deformation, which is a combination of the complete rotation matrix and the complete stretch tensor. It is also a local transformation. After image registration between the patient image and the atlas image, the deformation gradient or complete deformation gradient jF is computed. If jM represents the computed tensor from the femur atlas HRpQCT image, then tensor mapping based on Complete Deformation, termed here CD, is
    jMjFjMjF−1.
    (21)

Fabric atlas

In this section we present the methodology employed to construct a population-based atlas of fabric information. Differently from the diverse femur atlases described in section, a fabric atlas refers to a femur atlas model with a mean fabric tensor distribution. The overview of the construction of fabric atlas is presented in Fig 3. We follow a similar strategy as in cardiac DTI imaging for statistical analysis of cardiac fibres [21, 22]. Initially, a femur atlas HRpQCT image is chosen from the population. Image registration is performed between the femur atlas HRpQCT image IF and another femur HRpQCT image of the population IM. For each control point j of the femur atlas, the corresponding control point IM(T−1(j x)) is found, and fabric tensor jM is computed following the procedure described in previous section. The computed fabric tensors from all control points are then mapped to the femur atlas HRpQCT image. Mapping is performed by CR tensor mapping method, as it yielded best results compared to other tensor mapping methods (see result section). The same procedure is repeated for the rest of the femur HRpQCT images IM1, IM2, …IMN of the population and the respective fabric tensors jM1, jM2, ….jMN are computed. Then, the mean fabric tensor at each control point j is computed by averaging

jM¯=1Nn=1N(jRjMjRT).
(22)

The resulting mean fabric tensor is an arithmetic synthetic fabric tensor distribution that is mapped on the femur atlas HRpQCT image (MTP) being closest to the synthetic average femur, as described in previous section. Along with the other femur atlases presented in previous section, the resulting fabric atlas will be evaluated for prediction of patient femur fabric information, using the evaluation metrics presented in the next section.

Fig 3
Overview of the construction of fabric atlas (femur atlas model with mean fabric) is presented.

Materials and experiments

Datasource

The study was performed on a database of pairs of QCT and HRpQCT images of human proximal femora. The database comprises 36 pairs (17 males, 19 females with age 76±12 years, range 46–96 years) and were obtained from a previous study [4]. In summary, each femur was scanned with a calibration phantom (BDC Phantom, QMR Gmbh, Germany) in a clinical QCT (Brillance64, Phillips, Germany, intensity: 100 mA, voltage: 120 kV, voxel size: 0.33 × 0.33 × 1.00 mm3), and HRpQCT (Xtreme CT, Scanco, Switzerland, intensity: 900 μA, voltage: 60 kVp, voxel size: 0.082 × 0.082 × 0.082 mm3). The QCT images were rescaled to an isotropic voxel spacing (0.33 × 0.33 × 0.33 mm3) and were rigidly registered to the corresponding HRpQCT images. From the HRpQCT images, the cortical bone was masked out according to the procedure reported in [23].

Femur morphology

In order to assess how representative the selected database is with respect to the shape variability of the femur anatomy, a femur morphology study was first performed. To this end, an implicit coordinate system of the femur was constructed as shown in Fig 4. First, the femoral head center is defined by a mass center of a spherical region with maximal cross-section area. The neck axis is then computed by following the procedure reported by Kang et al. [24, 25]. In short, the radius of the spherical region of the femoral head is enlarged by 25%, and an initial neck center is defined. Using Powell’s optimization [26], the femoral neck center is computed, and the neck axis is defined as the line between femoral head center and femoral neck center (see Fig 4). The intersection point between the neck axis and the lateral surface of the femur is defined as the neck-axis-end-point. Then, the mass center of slices distal to this point are computed, followed by RANSAC fitting [27] to define the shaft axis. Generally, as the neck and shaft axes do not intersect, a mid point is defined as the shortest distance between the neck and shaft axes. The most distal point of the shaft axis is chosen as shaft-axis-distal-point. An implicit coordinate system is constructed by connecting femoral head center, mid point and shaft-axis-distal-point. As morphological parameters we calculated known shape descriptors of the femur, such as the caput-collum-diaphyseal angle (CCD), femoral head diameter, and distances describing the femoral neck anatomy. Femur morphology was computed for the total and sex-specific populations, and are summarized in Table 1. Between the Sex-specific populations, the morphology of the femurs were not found to be statistically significant(p>0.05).

Fig 4
Implicit coordinate system of the human femur.
Table 1
Morphology of the femur population.

Image pre-processing

Image pre-processing was performed on femur images to correct its shaft length, as the acquired images have varying shaft length. This step was also performed to ensure that the image registration step is not affected by differences in the anatomy. The shaft region of the femur was chopped such that the ratio between the distance femoral head center and mid point, and mid point and shaft-axis-distal-point (see Fig 4) equals 0.7, which was found empirically in order to yield an stable image registration. All the femurs were rigidly aligned with mid point as center.

Experimental design

We designed two experiments to answer the three open questions in registration-based bone fabric prediction, summarized below:

  1. Impact of femur atlas selection and sex considerations.
  2. Impact of the fabric mapping transformation on fabric prediction accuracy
  3. Potential of population-wide and sex-specific mean fabric atlases

In the first experiment, we combined in the evaluation the analysis of using different femur atlases (section) as well as different fabric tensor mapping transformations (section). In the second experiment, we evaluated the accuracy of predicting femur fabric by means of a femur atlas featuring a synthetically generated mean fabric (section) or its corresponding fabric tensor, as extracted from its HRpQCT image pair.

Evaluation scheme and metrics

For numerical evaluation a leave-one-out (LOO) strategy was followed. Specifically, a femur is chosen from the population as the patient’s femur and its counterpart (left or right) is removed from the population to remove bias in the analysis. At each control point j, the predicted fabric tensor for the patient’s femur QCT image is represented as M^ with eigenvalues m3^, m2^, m1^ and eigenvectors m3^, m2^, m1^ and the corresponding ground truth fabric tensor is computed from the patient’s femur HRpQCT image and is represented as M with eigenvalues m3, m2, m1 and eigenvectors m3, m2, m1. The process is repeated for each image in the database (N = 72 images), for each femur atlas (N = 6, section), and for each fabric mapping transformation (N = 5, section), leading to 5′040 (72 × 70) image registrations.

To evaluate the accuracy of the predicted femur fabric, we adopted the same evaluation metric as described in [6]. Namely, tensor norm error (TNerror), degree of anisotropy error (DAerror) and angular error of the principal tensor direction (PTDerror), are computed as follows:

TNerror=M^MM,DAerror=|DA^DA|DAandPTDerror=arccos(m3^,m3),
(23)

where the predicted, and ground-truth degree of anisotropy (DA) and (DA^), respectively, are computed as

DA^=m3^m1^andDA=m3m1.
(24)

The average error for each evaluation metric was computed for all the control points J and for all images in the LOO study, and were used as the base for comparison.

Results

Impact of femur atlas selection and sex considerations

We first present in Fig 5 overall results for all three evaluation metrics, for each femur atlas and fabric tensor mapping transformation. Regarding the selection of the femur atlas, as expected the farthest femur atlases (FTP and FSP) yielded the highest errors, followed by the mean femur atlases (MTP and MSP). We remark that the selection of FTP and FSP was motivated to reflect a potential worst-case scenario and to test the hypothesis that an atlas should be as similar as possible to the patient image on which fabric is predicted. Results on all metrics showed that choosing the closest femur atlases (CTP, CSP) yields the lowest errors, which verifies the importance of the femur atlas selection.

Fig 5
Summary of prediction error for different sets of combination of tensor mapping methods and femur atlases.

Regarding sex, no statistically significant differences for all three metrics were found (p>0.05) between choosing atlases from the total population (CTP, FTP, MTP) or sex-specific ones (CSP, FSP, MSP). This results suggests that it might not be necessary to create sex-specific femur atlases when predicting femur fabric.

Impact of the fabric mapping transformation on fabric prediction accuracy

Regarding the impact of the fabric mapping transformation, results presented in Fig 5 show that fabric tensor mapping methods involving only rotation components (CR,AR) produce lower errors than tensor mapping methods involving both rotation and stretch components (CD,AD). Among the methods relying only on rotation, fabric tensor mapping by CR yielded the lowest error, followed by AR and NR. However, only the TNerror and PTDerror were found to be significantly different, as shown in Fig 6. This is due to the fact that fabric tensor mapping methods involving only a rotation component do not alter the eigenvalues, and hence DAerror remained the same.

Fig 6
Summary of prediction error for different sets of combination of fabric tensor mapping methods involving only rotation component (NR, AR, CR) and femur atlases(CTP in blue, MTP in pink).

Relative to the selected femur atlas, using CR fabric tensor mapping in combination with CTP yielded TNerror = 7.3±0.9%, DAerror = 6.6±1.3%, and PTDerror = 25±2°). These results compare favorably to those yielded when using MTP as femur atlas, with TNerror = 7.7±1.0%, DAerror = 7.0±1.4%, and PTDerror = 25±2°. Nonetheless, it is remarked that while CTP requires image registration for each image of the database to calculate the distance metric DM, MTP is computed only once and does not require further computations across the database. Fig 6 focuses on analyzing the performance for these two femur atlases (CTP,MTP), and fabric tensor mapping methods based on a rotation component (NR, AR, CR). Statistically significant differences (p>0.05) were found between NR and CR, and between CR and AR, but not between NR and AR, confirming the value of using CR as preferred fabric tensor mapping transformation.

Potential of population-wide and sex-specific mean fabric atlases

Fig 7, shows fabric prediction errors for all three evaluation metrics, when predicting femur fabric by means of a femur atlas featuring a synthetically generated mean fabric or by its corresponding real fabric tensor, as extracted from its HRpQCT image pair. In this experiment, MTP and CR were chosen as femur atlas and fabric tensor mapping method, respectively. We found that using the synthetically generated fabric atlas yielded higher error than using the real fabric from the corresponding HRpQCT fabric. A statistical difference was found (p < 0.05).

Fig 7
Summary of prediction error, comparing bone fabric prediction by femur atlas with HRpQCT fabric and fabric atlas.

Spatial and bone mineral density based evaluation of femur fabric prediction

We performed a spatial analysis of fabric prediction performance to analyze how the prediction errors are spatially distributed. Fig 8, shows in three different planes, the prediction of bone fabric for an example case using selected femur atlases CTP and MTP, and CR as fabric tensor mapping method. It is observed that the TNerror varies widely across different regions of the femur. We observed that lower error are observed across the main loading direction and in femoral head regions. Higher errors were observed in the shaft and in lower trochanter regions.

Fig 8
Illustration of bone fabric prediction accuracy achieved on a test case femur.

Finally, as the registration process is driven by image intensity information we were interested to analyze whether there is a correlation between bone mineral density and fabric prediction error. Fig 9 shows for each metric bone fabric prediction errors for different Bone Volume over Total Volume (BVTV) bins. In this experiment, MTP and CR were chosen as femur atlas and fabric tensor mapping method, respectively. We observed increasing errors for TNerror and DAerror in regions of moderate to high BVTV, whereas lower PTDerror errors were found for moderate to high BVTV regions.

Fig 9
Summary of prediction error comparing bone fabric prediction for different BVTV bins.

Discussion

In this study we propose a novel registration-based estimation of bone fabric directly from clinical QCT images. It is designed to preserve tensor properties of bone fabric and to map bone fabric by a global and local decomposition of the gradient of a non-rigid image registration transformation. We analyzed the importance of the fabric tensor mapping transformation as well as the femur atlas used to map the fabric information into a target QCT image. We further evaluated and reported the performance of a population-, and sex-specific atlas of bone fabric used within the proposed registration-based fabric estimation approach.

The entire study was performed on a database of 36 pairs of human proximal femora [4], for which the results of the morphology analysis suggest that the femurs used in the present study are representative of femurs from other studies [28].

Importance of femur atlas selection

Regarding different femur atlases, it becomes clear from Fig 5 that the farthest to the patient’s femur, FTP and FSP, yielded rather poor results. Conversely, the closest to the patient’s femur, CTP and CSP, yielded the best results, which allow us to conclude that bone fabric prediction based on image registration is sensitive to the selected femur atlas. These results are in agreement with the strategy presented in [3] where a femur database and selection scheme was originally presented. From a physiological loading point of view, it is indeed expected that differences in bone anatomy have an effect on the underlying bone fabric [29]. As reported in Table 1 as well as in previous studies regarding bone femur morphology [28], such difference in bone anatomy is observed through parameters such as the CCD angle and neck length. However, further FE simulations on a representative population are required to assess the impact of femur atlas selection on bone strength prediction.

In addition, regarding sex considerations, results suggests that there is no major benefit in employing sex-specific femur atlases for fabric prediction. The probable reason for this finding is that the variability in femoral shape (in terms of DM) between sex-specific populations is smaller than 3% of the population shape variability. Table 1 supports this statement where the femurs’ morphological variables of females and males are in average similar (p>0.05).

Interestingly, results presented in Figs Figs55 and and66 suggest while the highest fabric prediction accuracy is attained with CTP, followed by MTP, their differences in accuracy are often statistically significant, but quantitatively the results are rather close. In this regard, one important practical limitation of using CTP involves computing the closest femur image (in terms of DM) in the population. In practice, such computations are prohibitive for large databases. On the contrary, MTP is computed once and if needed, it can be updated for an extended or different population database.

Impact of fabric tensor mapping

Looking at different fabric tensor mapping methods, results indicate that tensor mapping by AR and CR performs better than CD and AD. The DA error clearly supports the conclusion that the stretch component of CD and AD tends to alter the bone fabric excessively. The importance of the tensor mapping method becomes also clear from Fig 6, where the tensor mapping by NR does not improve fabric prediction in terms of TN and PTD (p<0.05). Conversely, fabric tensor mapping by AR does improve fabric prediction (p<0.05), and fabric tensor mapping by CR yields the best fabric prediction accuracy (p<0.05).

On the other hand, as rotation-only mapping approaches do not alter the eigenvalues of fabric tensors, these approaches are not capable of predicting DA. These results suggest that the degrees of freedom of the chosen transformation model plays an important role, and a trade-off between accuracy of predicting fabric orientation and DA needs to be considered when using rotation-only mapping schemes.

Turning to the concept of employing a synthetically generated fabric atlas, results presented in Fig 7 suggests that bone fabric predictions are considerably less accurate than when using the real fabric of the HRpQCT image. One possible reason of this is the fact that the femur atlas stems from a real bone image, which is then combined with a computed (i.e synthetic) mean fabric distribution. Such combination might not fully characterize the interplay between bone morphology and fabric distribution, as naturally occurs for a real bone image, where bone fabric and bone morphology are interrelated [29]. Although algorithms exist to compute a mean femur atlas [30], our experiments yielded an over-smoothed synthetic image not preserving the required image quality inherent of HRpQCT. Further research on atlas construction approaches specifically designed to deal with tensorial information, such as [21, 22], might provide improvements to the creation of high resolution femur atlases.

Spatial and BVTV based analysis of fabric prediction accuracy

The spatial distribution of the tensor norm error in Fig 8 shows that the bone fabric prediction accuracy varies widely across regions. In particular, the femoral head and the main loading trajectory (which is of primary interest for FE analysis) present higher fabric prediction accuracy than the inter-trochanteric region. This finding may be due to the inability of the image registration approach to handle properly regions of lower BV/TV. From a physiological point of view, the higher BV/TV relates to bone micro-architecture oriented along the principal stresses acting on the femur and forming characteristic trajectories [31]. To further clarify this issue, the bone fabric prediction was analyzed with respect to BV/TV. Fig 9 indicates indeed better predictions of the major fabric orientation in regions of high BV/TV.

Limitations of the present study

Some limitations of this study have to be mentioned. First, we use a distance metric (DM) for selection of the femur atlas that does not not consider any anthropometric parameters or ethnic variation. Their inclusion in the analysis would be of great interest for patient-specific FE analysis. However, such information was not available for the database used in this study.

Second, the processing time of the present approach for one femur takes 40 min. The time was measured on a desktop with the application running single-threaded on a 3.20 GHz Intel Core i7 processor. Such computation time is relatively high compared to other machine learning based approaches for bone fabric predictions [6, 7]. However, this is a major common disadvantage of all image-based registration approaches.

Third, the optimal mapping approach, CR, is not capable of improving the prediction of DA. One of the possible way to address this problem would be the use of poly-affine registration [32, 33], where a set of affine transformations is employed to characterize spatial transformations with a low number of parameters. As described, our experiments suggest that a trade-off between flexibility of the transformation model to morph the atlas image onto the patient image, and preservation of bone fabric information exists. Similarly, the use of dedicated registration algorithms encoding specific properties linked to the anatomy or disease in study (e.g. [34]), or registration approaches previously proposed for Diffusion Tensor Imaging (DTI) (e.g. [21, 22]) might provide a better fabric prediction based on image-registration approaches.

Comparison to previous approaches

Although previous approaches have used different datasets and sample sizes (Taghizadeh et.al [2] N = 10, Chandran et.al [7] N = 30, Lekadir et.al [6] N = 33), some quantitative comparisons are worth mentioning. The present study yielded lower TN and PTD errors than the one of [2], where a TN error of 14.8 ± 1.5% and a PTD error of 29.7 ± 3.3° were reported. On the other hand, studies based on machine learning approaches, such as [7] reported a TN error of 6 ± 2%, and a PTD error of 19 ± 7°, and [6] reported a TN error of 7 ± 1% and a PTD error of 15.6 ± 2.3°, which are comparable for TN but lower for PTD compared to the registration-based method explored in the present study. Regarding prediction of DA, the studies of [7] and [6] reported DA error (6 ± 2% and of 7 ± 1% respectively), which are comparable to the ones obtained with registration-based methods.

The present study lacks the experimental data to validate the role of predicted bone fabric in computational models for calculation of bone strength. However, previously reported bone fabric prediction accuracies [3, 6, 29, 35] are in the similar range of prediction accuracy reported here. Hence, we expect corresponding improvements in bone strength predictions.

Conclusion

In conclusion, we proposed a novel image-registration based femur fabric prediction directly from clinical QCT image. The methodology is robust and favorably compares to previous state of the art registration-based method for femur fabric prediction. Furthermore, we present a comprehensive analysis of key components of the registration-based approach for bone fabric prediction in the proximal femur. From the results, we could answer three open questions. First, compromising between accuracy and computing time, the optimal femur atlas corresponds to the mean of the total population (MTP). Second, the best tensor mapping method is provided by complete rotation (CR). Third, a population average fabric atlas produced higher errors in fabric prediction than employing directly MTP and CR, and hence it is not recommended. By employing MTP, registration with a whole database of femurs becomes unnecessary and reduces considerably computational time.

The reported findings are promising for a clinical implementation and exploitation for patient-specific analysis as it is has potential to leverage bone architectural information directly from standard clinical imaging. Moreover, while image registration algorithms are improving we note on the importance of designing clinically- and task-oriented image registration pipelines. In this sense, the set of recommendations generated from this study are expected to guide the development of dedicated image based assessment methodologies of bone architecture from clinical imaging. The impact of the identified image-registration methodology on the prediction of hip strength by finite element analysis will be evaluated in future work.

Acknowledgments

This work is supported by the Swiss National Science Foundation, Project number—143769. The authors thank Dr. Ghislain Maquer for his comments on mechanical aspects and proof reading. The authors would like to thank Enrico Dall’Ara and Dieter H. Pahr (TU Wien, Austria) for sharing their QCT and HRpQCT images.

Funding Statement

This work is supported by the Swiss National Science Foundation (http://p3.snf.ch/Project-143769.

Data Availability

Data Availability

Data used in our work are third-party data obtained from the following source: Dall’Ara E, Luisier B, Schmidt R, Kainberger F, Zysset P, Pahr D. A nonlinear QCT-based finite element model validation study for the human femur tested in two configurations in vitro. Bone. 2013 Jan;52(1):27-38. The data will be available upon request to all interested researchers. Please contact PD DI. Dr. Dieter H. PAHR (Institute of Lightweight Design and Structural Biomechanics (E317)) at ta.ca.neiwut.bsli@rhap to request data access.

References

1. Taghizadeh E, Chandran V, Reyes M, Zysset P, Büchler P. Statistical analysis of the inter-individual variations of the bone shape, volume fraction and fabric and their correlations in the proximal femur. Bone. 2017. October;103:252–261. [PubMed]
2. Taghizadeh E, Reyes M, Zysset P, Latypova A, Terrier A, Büchler P. Biomechanical Role of Bone Anisotropy Estimated on Clinical CT Scans by Image Registration. Ann Biomed Eng. 2016. August;44(8):2505–17. doi: 10.1007/s10439-016-1551-4 [PubMed]
3. Hazrati-Marangalou J, Ito K, Cataldi M, Taddei F and van Rietbergen B. A novel approach to estimate trabecular bone anisotropy using a database approach. J Biomech. 2013. September 27;46(14):2356–62. doi: 10.1016/j.jbiomech.2013.07.042 [PubMed]
4. Dall’Ara E, Luisier B, Schmidt R, Kainberger F, Zysset P, Pahr D. A nonlinear QCT-based finite element model validation study for the human femur tested in two configurations in vitro. Bone. 2013. January;52(1):27–38. doi: 10.1016/j.bone.2012.09.006 [PubMed]
5. Luisier B, Dall’Ara E, Pahr DH. Orthotropic HR-pQCT-based FE models improve strength predictions for stance but not for side-way fall loading compared to isotropic QCT-based FE models of human femurs. J Mech Behav Biomed Mater. 2014. April;32:287–299. doi: 10.1016/j.jmbbm.2014.01.006 [PubMed]
6. Lekadir K, Hazrati-Marangalou J, Hoogendoorn C, Taylor Z, van Rietbergen B, Frangi AF. Statistical estimation of femur micro-architecture using optimal shape and density predictors. J Biomech. 2015. February 26;48(4):598–603. doi: 10.1016/j.jbiomech.2015.01.002 [PubMed]
7. Chandran V, Zysset P, Reyes M. Prediction of Trabecular Bone Anisotropy from Quantitative Computed Tomography Using Supervised Learning and a Novel Morphometric Feature Descriptor. Med Image Comput Comput Assist Interv. 2015. November; 9349:621–628.
8. Taghizadeh E, Maquer G, Reyes M, Büchler P. Including the trabecular anisotropy from registered microCT data in homogenized FE model improves the bone’s mechanical predictions. CMBBE. 2014.
9. Bonaretti S, Seiler C, Boichon C, Reyes M, Büchler P. Image-based vs mesh-based statistical appearance models of the human femur: Implications for finite element simulations. Med Eng Phys. 2014. December;36(12):1626–35. doi: 10.1016/j.medengphy.2014.09.006 [PubMed]
10. Klein S, Staring M, Murphy K, Viergever MA, Pluim JP. elastix: a toolbox for intensity-based medical image registration. IEEE Trans Med Imaging. 2010. January;29(1):196–205. doi: 10.1109/TMI.2009.2035616 [PubMed]
11. Rueckert D, Sonoda LI, Hayes C, Hill DL, Leach MO, Hawkes DJ. Nonrigid registration using free-form deformations: application to breast MR images. IEEE Trans Med Imaging. 1999. August;18(8):712–21. doi: 10.1109/42.796284 [PubMed]
12. Dice LR. Measures of the Amount of Ecologic Association Between Species. Ecology. 1945;26(3):297–302. doi: 10.2307/1932409
13. Elhawary H, Oguro S, Tuncali K, Morrison PR, Tatli S, Shyn PB, Silverman SG, Hata N. Multimodality non-rigid image registration for planning, targeting and monitoring during CT-guided percutaneous liver tumor cryoablation. Acad Radiol. 2010. November;17(11):1334–44. doi: 10.1016/j.acra.2010.06.004 [PMC free article] [PubMed]
14. Fréchet M. Les éléments aléatoires de nature quelconque dans un espace distancié. Annales de l’institut Henri Poincaré. 1948;10(4):215–310.
15. Lepore N, Brun C, Chou YY, Chiang MC, Dutton RA, Hayashi KM, Luders E, Lopez OL, Aizenstein HJ, Toga AW, Becker JT, Thompson PM. Generalized tensor-based morphometry of HIV/AIDS using multivariate statistics on deformation tensors. IEEE Trans Med Imaging. 2008. January;27(1):129–41. doi: 10.1109/TMI.2007.906091 [PMC free article] [PubMed]
16. Joshi S, Davis B, Jomier M, Gerig G. Unbiased diffeomorphic atlas construction for computational anatomy. Neuroimage. 2004;23 Suppl 1:S151–60. doi: 10.1016/j.neuroimage.2004.07.068 [PubMed]
17. Kalinic H. Atlas-based image segmentation: A Survey. Croatian Scientific Bibliography. 2009;1–7.
18. Harrigan TP, Mann RW. Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor. Journal of Materials Science. 1984;19(3):761–767. doi: 10.1007/BF00540446
19. Varga P, Zysset P. Assessment of volume fraction and fabric in the distal radius using HR-pQCT. Bone. 2009. November;45(5):909–17. doi: 10.1016/j.bone.2009.07.001 [PubMed]
20. Geng X, Styner M, Gupta A, Shen D, Gilmore JH. Multi-contrast diffusion tensor image registration with structural MRI. ISBI. 2012. May;684–687.
21. Peyrat JM, Sermesant M, Pennec X, Delingette H, Xu C, McVeigh ER, Ayache N. Towards a statistical atlas of cardiac fiber structure. Med Image Comput Comput Assist Interv. 2006;9(Pt 1):297–304. [PubMed]
22. Peyrat JM, Sermesant M, Pennec X, Delingette H, Xu C, McVeigh ER, Ayache N. A computational framework for the statistical analysis of cardiac diffusion tensors: application to a small database of canine hearts. IEEE Trans Med Imaging. 2007. November;26(11):1500–14. doi: 10.1109/TMI.2007.907286 [PubMed]
23. Pahr DH, Zysset P. From high-resolution CT data to finite element models: development of an integrated modular framework. Comput Methods Biomech Biomed Engin. 2009. February;12(1):45–57. doi: 10.1080/10255840802144105 [PubMed]
24. Kang Y, Engelke K, Fuchs C, Kalender WA. An anatomic coordinate system of the femoral neck for highly reproducible BMD measurements using 3D QCT. Comput Med Imaging Graph. 2005. October;29(7):533–41. doi: 10.1016/j.compmedimag.2005.05.002 [PubMed]
25. Kang Y, Engelke K, Kalender WA. A new accurate and precise 3-D segmentation method for skeletal structures in volumetric CT data. IEEE Trans Med Imaging. 2003. May;22(5):586–98. doi: 10.1109/TMI.2003.812265 [PubMed]
26. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press; 2007;3rd ed.
27. Fischler MA, Bolles RC. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM. 1981. June;24(6):381–395. doi: 10.1145/358669.358692
28. Bergot C, Bousson V, Meunier A, Laval-Jeantet M, Laredo JD. Hip fracture risk and proximal femur geometry from DXA scans. Osteoporosis International. 2002;13(7):542–550. doi: 10.1007/s001980200071 [PubMed]
29. Hazrati Marangalou J, Ito K, van Rietbergen B. A novel approach to estimate trabecular bone anisotropy from stress tensors. Biomech Model Mechanobiol. 2015. January;14(1):39–48. doi: 10.1007/s10237-014-0584-6 [PubMed]
30. Rueckert D, Frangi AF, Schnabel JA. Automatic Construction of 3D Statistical Deformation Models Using Non-rigid Registration. Med Image Comput Comput Assist Interv. 2001;9(Pt 1):77–84.
31. Wolff J. Das Gesetz der Transformation der knochen. AHirchwild. 1892.
32. Seiler C, Pennec X, Reyes M. Geometry-aware multiscale image registration via OBB Tree-based polyaffine log-demons. Med Image Comput Comput Assist Interv. 2011;14(Pt 2):631–638. [PubMed]
33. Seiler C, Pennec X, Reyes M. Simultaneous Multiscale Polyaffine Registration by Incorporating Deformation Statistics. Med Image Comput Comput Assist Interv. 2012;15(Pt 2):130–137. [PubMed]
34. McLeod K, Sermesant M, Beerbaum P, Pennec X. Spatio-Temporal Tensor Decomposition of a Polyaffine Motion Model for a Better Analysis of Pathological Left Ventricular Dynamics. IEEE Trans Med Imaging. 2015. February 19 doi: 10.1109/TMI.2015.2405579 [PubMed]
35. Lekadir K, Noble C, Hazrati-Marangalou J, Hoogendoorn C, van Rietbergen B, Taylor ZA, Frangi AF. Patient-Specific Biomechanical Modeling of Bone Strength Using Statistically-Derived Fabric Tensors. Ann Biomed Eng. 2016. January;44(1):234–46. doi: 10.1007/s10439-015-1432-2 [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science