High resolution cortical bone thickness measurement from clinical CT data

Abstract

The distribution of cortical bone in the proximal femur is believed to be a critical component in determining fracture resistance. Current CT technology is limited in its ability to measure cortical thickness, especially in the sub-millimetre range which lies within the point spread function of today's clinical scanners. In this paper, we present a novel technique that is capable of producing unbiased thickness estimates down to 0.3mm. The technique relies on a mathematical model of the anatomy and the imaging system, which is fitted to the data at a large number of sites around the proximal femur, producing around 17,000 independent thickness estimates per specimen. In a series of experiments on 16 cadaveric femurs, estimation errors were measured as -0.01+/-0.58mm (mean+/-1std.dev.) for cortical thicknesses in the range 0.3-4mm. This compares with 0.25+/-0.69mm for simple thresholding and 0.90+/-0.92mm for a variant of the 50% relative threshold method. In the clinically relevant sub-millimetre range, thresholding increasingly fails to detect the cortex at all, whereas the new technique continues to perform well. The many cortical thickness estimates can be displayed as a colour map painted onto the femoral surface. Computation of the surfaces and colour maps is largely automatic, requiring around 15min on a modest laptop computer.

Fig. 1

The effect of surface orientation and CT slice thickness on the apparent thickness and out-of-plane PSF. A finite CT slice width s causes cortical layers oblique to the slice to appear as if they have been blurred with a rectangular PSF, with extent dependent on surface orientation a, leading to a piecewise linear image. The true cortical thickness $t_r$ can be expressed as a function of a and the measured thickness $t_m$.

Fig. 2

A simulation of the cortical imaging process. The CT value along a line passing through the cortex is simulated by blurring a piecewise constant density function (with CT values of −1000, 1750 and 0) with a Gaussian kernel of standard deviation 1 mm. The solid lines show the underlying and blurred CT values, the dotted lines show a threshold of 600, and the dashed lines show values half way between the appropriate CT baseline and the blurred peak (50% relative threshold). In (a), the cortex is sufficiently wide for thresholding and the 50% relative threshold methods to work correctly. In (b), the threshold still gives approximately the correct thickness, but the 50% relative threshold method overestimates. In (c), no values exceed the threshold, and the 50% relative threshold method now reports the width of the PSF rather than that of the cortex.

Fig. 3

Cortical thickness in high and low resolution CT data. (a) and (b) show corresponding CT slices from a cadaveric femur scanned at high resolution (Xtreme pQCT, 82 μm/pixel, Scanco Medical AG, Brüttisellen, Switzerland) and low resolution (Siemens Somatom Sensation 64 MDCT, 589 μm/pixel, Siemens AG, Erlangen, Germany). (c) and (d) show the CT values along the lines labelled (1) and (2) in the images. In (c), the cortical layer is sufficiently thick to show a significant peak in the low resolution data. In (d), the thickness is somewhat less, and in this case there is only a very shallow peak in the low resolution data. Also shown in (c) and (d) are the best fit models $y_{\text{blur}}$, and the corresponding piecewise constant density functions.

Fig. 4

Discrimination of cortical and trabecular bone peaks for the half-max and threshold methods. (a) and (b) are two examples from the high resolution data showing how the model-fitting process establishes which peaks are within the cortex. The upper plots correspond to the half-max method, the lower plots to the threshold method. The model-fitting algorithm optimally (in a least squares sense) divides the data into spans of higher (cortical) and lower (trabecular, exterior) density values. For the half-max method, the boundaries tend to be positioned very close to the half way point between the cortical and trabecular/exterior values. For the threshold method, the boundaries are positioned precisely on one of the threshold crossing points.

Fig. 5

Creation of the surface and the use of surface normals for guiding thickness estimation. A surface is generated by thresholding the entire data set, then extracting contours in each plane to sub-pixel resolution. Contours are then locally edited (a) to correct erroneously excluded regions and (b) to remove any adjoining structures. A surface is interpolated through these contours, and the surface vertices and normals used to guide the in-plane thickness estimation. (c) The surface normals are shown in cyan, cortical edges as red dots, and the plot shows the thickness estimation process for the red line at the bottom.

Fig. 6

Mapping and filtering of high resolution thickness estimates. (a) Once cortical thickness has been estimated at each vertex, this is mapped back onto the surface as a colour. (b) The high resolution thickness map is filtered over the surface, such that the spatial resolution of the map is approximately equivalent to that of the low resolution data.

Fig. 7

Co-registration of high and low resolution surfaces. High resolution (a) and low resolution (b) surfaces are aligned using manual registration followed by the ICP algorithm. (c) The registered surfaces. (d) Point-wise alignment discrepancies are generally well below 1 mm.

Fig. 8

Distribution of cortical thickness and errors across all 16 cadaveric data sets. The top graph shows the thickness distribution derived from the high resolution data. The bottom graph shows the error distribution for the ‘new variable’ measurement method. Despite the narrow peak, there are outlying errors up to ±15 mm. The dotted lines indicate the 4 mm limits used to exclude poor matches (less than 3% of the total) from some of the subsequent error analysis.

Fig. 9

Difficulties comparing low and high resolution thickness measurements. (a) and (b) show corresponding high and low resolution slices through the CT data. (c) and (d) are plots of the CT values along the lines labelled (1) and (2) respectively. In (c), the surface normal in the high resolution data is not compatible with that of the low resolution data. In (d), there is ambiguity in the definition of cortical thickness: both estimators detect the outermost resolvable layer, but this layer is not the same at the two resolutions.

Fig. 10

Thickness estimation on simulated data with different values of cortical thickness in the range 0.2–5.5 mm. The CT value along a line passing through the cortical layer was simulated by blurring a piecewise constant density function with a Gaussian kernel of standard deviation 1.25 mm. The unblurred cortical density was set to 1750 HU and the trabecular density to 0 HU. The cortical thickness was then estimated using a variety of techniques: the half-max method, thresholding at different CT values from 300 HU to 1700 HU, and the new method assuming cortical CT values $(y_1)$ in the range 1300 HU to 2200 HU. The results are presented as follows. There is a single dot-dashed line for the half-max method, one dotted line for each threshold value (thresholding method), and one dashed line for each assumed cortical CT value (new method). (a) shows a simulated scan in air (unblurred extra-cortical density of −1000 HU) and (b) shows an in vivo simulation (unblurred extra-cortical density of 0 HU). In both cases, the half-max method overestimates at low thicknesses. Thresholding gives better results for carefully chosen thresholds: note how the optimal threshold varies significantly from 600 HU in (a) to 800 HU in (b). The new method gives consistently good results, with relatively little sensitivity to the assumed cortical CT value: the dashed lines are grouped closely around the main diagonal.

Fig. 11

Cadaveric cortical thickness estimation as a function of thickness. The upper graph in each case shows the mean ± 1 standard deviation of $t_l$ on the y-axis, for each 0.1 mm band of $t_h$ on the x-axis. This is similar to the presentation used for the simulated data in Fig. 10. The lower graph indicates the percentage of successful thickness estimates in each band. (a) shows results for the best case (473), (b) for the worst (467) and (c) the combined results across all 16 femurs. The three colours have been rendered with some transparency, so the ±1 standard deviation ranges remain distinguishable where they overlap. At low thicknesses, the half-max and threshold results are not what one would expect from the model (Fig. 10). We discuss possible reasons for this in Section 4.

Fig. 12

Thickness estimation for the best case (473). (a)–(d) are thickness maps from the front, (h)–(k) from the back. (e)–(g) and (l)–(n) are absolute thickness error maps. The first column (a, h) shows the high resolution ‘ground truth’ measurements. The second column (b, e, i, l) shows the thresholding results, the third column (c, f, j, m) half-max, and the fourth column (d, g, k, n) the ‘new variable’ results.

Fig. 13

Thickness estimation for the worst case (467). (a)–(d) are thickness maps from the front, (h)–(k) from the back. (e)–(g) and (l)–(n) are absolute thickness error maps. The first column (a, h) shows the high resolution ‘ground truth’ measurements. The second column (b, e, i, l) shows the thresholding results, the third column (c, f, j, m) half-max, and the fourth column (d, g, k, n) the ‘new variable’ results.

Fig. 14

Cortical thickness estimation for an 87-year old male volunteer with a sub-capital femoral neck fracture. The maps were produced using the ‘new variable’ method. (a) and (b) show the contralateral and fractured femurs from the back, (c) and (d) from the front. This demonstrates the feasibility of estimating cortical thickness maps in vivo, even for fractured femurs. The fractured and contralateral maps are reassuringly similar.

Fig. 15

Cortical thickness estimation at 1 mm and 3 mm CT slice thickness for a 39-year old female volunteer. The maps were produced using the ‘new variable’ method. (a) and (b) are estimates from CT data reconstructed at 1 mm slice thickness, (c) and (d) from the same data reconstructed at 3 mm slice thickness. (e) and (f) are absolute difference maps, revealing very few discrepancies between the two sets of results.