A novel interpolation approach for the generation of 3D-geometric digital bone models from image stacks.

The authors propose a new 3D interpolation algorithm for the generation of digital geometric 3D-models of bones from existing image stacks obtained by peripheral Quantitative Computed Tomography (pQCT) or Magnetic Resonance Imaging (MRI). The technique is based on the interpolation of radial gray value profiles of the pQCT cross sections. The method has been validated by using an ex-vivo human tibia and by comparing interpolated pQCT images with images from scans taken at the same position. A diversity index of ⟨0.4 (1 meaning maximal diversity) even for the structurally complex region of the epiphysis, along with the good agreement of mineral-density-weighted cross-sectional moment of inertia (CSMI), demonstrate the high quality of our interpolation approach. Thus the authors demonstrate that this interpolation scheme can substantially improve the generation of 3D models from sparse scan sets, not only with respect to the outer shape but also with respect to the internal gray-value derived material property distribution.

Introduction

Reasonable simulations with Finite Element (FE) techniques require useful 3D-models. Such customized 3D-models of long bones like the tibia or femur can be calculated from cross-sectional images obtained in vivo by peripheral Quantitative Computer Tomography (pQCT).

Single-slice pQCT machines (e.g. Stratec XCT series) are used in many labs around the world. These machines are small, involve low radiation doses and yield highly valuable information about volumetric bone mineral density. However, the acquisition time for one cross-sectional image is typically about 90 seconds. Therefore it is not easy to scan an entire tibia bone with slice interspacing of less than 1 cm. In consequence, data sets are typically obtained with larger inter-slice distances. This makes the computation of a reasonable 3D bone model with adequate interslice 2D resolution difficult. Ideally, a 3D-model consists of isometric voxels at an internal resolution of about 0.5 mm typical of the resolution of conventional pQCT.

In an earlier study we had worked with slice distances of 1-2 cm and were looking for adequate interpolation mechanisms in order to apply FE simulation to further exploit the data of that study. Commercial software packages (e.g. MIMICS Materialise, Belgium) provide tools for generating 3D-models out of a series of cross-sectional scans. This kind of software is able to generate 3D surface models based on the surface pattern under a given gray-value threshold. It also allows for smoothing of the model surfaces and thereby interpolation of the outer shape. However, the gray-value distribution within the surface model becomes distorted during this process or even gets lost, as the contour shaping functions ignore the gray-values inside. Moreover, the quality of the shape smoothing depends very strongly on the shape transition from slice to slice: while the smoothing might be acceptable for the steady changes in the mid-shaft of a long bone it produces a terraced design for the curved surfaces of the epiphyses.

Typical stack of pQCT-scans of the proximal part of a human tibia; slice distance 10 mm.

terraced shape of a surface model generated with MIMICS (MIMICS Research v 18.0, Materialize) from 10 mm spaced pQCT scans of the proximal part of a human tibia; A: isometric view, B: a coronal sectional view.

Therefore, we looked into existing interpolation methods. The work on image processing problems has produced a wealth of interpolation techniques[1], but most of them deal with interpolation of features inside an image. There are much fewer approaches dealing with 3D inter-image interpolation. Thus, shape based interpolation based on surface evolution, but neglecting the gray value distribution has been used for bone shape interpolation[2], image based dosimetric models[3], organ and tissue discrimination[4,5], tissue feature displacement[6] and CT scanning planning[7]. Other techniques use classical interpolation algorithms, such as the ‘nearest neighbour’ approach, trilinear or B-spline techniques in combination with segmentation to a bone-no bone threshold, applied for example for the quantitative bone morphometry in the context of the registration of time-lapsed micro-CT scans[8]. All those approaches focus more or less on the surface shape and neglect any gray-value distribution inside the volume of the given structure. We found in preliminary work that classical interpolation alone, without doing any bone to no-bone segmentation, results in acceptable interpolates for the gray-value distribution inside the bone volumes but fails completely when it comes to a shift of sharp bone edges. Initial tests with the Matlab (The MathWorks Inc., Natick, MA) function “interpn” for multidimensional interpolation showed that the result introduces discontinuities of the bone surface that are not physically meaningful. Even in a more recent work that used an advanced projection-on-convex-curve-set (POCS) technique[9], this principle problem seems unsolved when looking at the presented interpolates (of MRI images in this case) which show a pronounced terrace pattern at the tissue boundaries. Overlaying a later segmentation against a threshold (as in[7]) might mitigate that behaviour sufficiently, but then, as a downside, the gray-value distribution is lost. We also looked into morphing techniques[10] but decided that the available methods did not really address the special problems of our case. Thus we had to find a solution that combines a shape interpolation with an interpolation of gray-values. The solution should be simple and easy to implement and not as complex as a Kriging algorithm[11]. In this paper we present a novel algorithm that takes radial profiles as the baseline for interpolation.

Methods

As stated before, the principle of the method is to perform two interpolations at the same time: an interpolation of the outer shape of a bone as well as an interpolation of the internal gray-value distribution representing the inner structure of the bone and its local mineral density. This is achieved by taking the radial profile of this entity as the baseline element for interpolation.

Figure 3 depicts the principles of the method. The term ‘radial profile’ means the (one-dimensional) distribution of gray values along a radial ray starting at a point near the center of a cross-section of a circular bone structure. The principle idea is quite straightforward: A complete set of radial profiles for the entire 360° of a given cross section scan of a long bone provides us with a good description of:

Figure 3

Scheme of the linear interpolation; the case for an interpolation factor of 2 (only one interpolated slice at half distance) is depicted; red cross: manually selected center; yellow cross: manually selected radius of the interpolation area; left side: first base slice with one of the radial rays in blue; right side second base slice; corresponding radial ray in green. The dimension of the radial rays is distance from center in number of voxels. Depending on the direction of the ray w.r.t. the square shaped voxels the absolute distance is varying. The arrows show the shifting of the radial rays during the procedure (see in the text).

the outer shape by the radial position of the onset of bone density (a distinct jump in density when approaching the bone from the outside) as well as of

the internal structure by the profile of gray-values along the radial rays.

That allows for interpolation of the shape by interpolating the position of the bone onset and the internal structure by interpolating the gray-value-profile taking the bone onset point as a reference point. Performing this for a complete set of radial rays an interpolated slice is built up.

The interpolation technique itself can utilize linear or cubic spline interpolation. An arbitrary factor can be selected that defines how many equidistant slices should be generated between the two base slices.

[Figure 3] shows the situation for the more simple linear interpolation approach. In this case the interpolation procedure is as follows:

Select two adjacent slices as base slices for the interpolation:

Provide the original distance between the base slices and the desired distance between the interpolated slices; an interpolation factor will be calculated; e.g. if it is needed to have interpolated slices with a distance of 2 mm starting from 10 mm distant base slices the interpolation factor is 5.

Manually: Select a center point in one of the base slices (red cross in Figure 3); it is not necessary to hit the real center correctly.

Manually: Select a point outside the bone cross-section (e.g. yellow cross in Figure 3); the distance between the center point and this second point will define the radius of a circle around the center point; take care that the entire bone cross-section of both base slices is inside this circle. It defines the maximum length of a radial ray and the interpolation area; all values outside this circle will not be considered.

Take one of the base slices as the baseline of the new interpolated slice. Inside a 360 degrees loop with an arbitrary angle step size (normally 1° or 0.5°) around the center point, and for every angle (and corresponding radial ray) do the following steps:

Generate the two radial profiles for the base slices (blue for slice #1, green for slice #2 in Figure 3).

Define the onset point of bone related gray-value (the outer edge of the bone; see explanation in Figure 3). These points will be reference points for the following interpolation processes.

Perform linear interpolation of bone onset: the bone onset for the interpolated profile is linearly interpolated between the two reference points for the new target slice taking interpolation factor and targeted position into account (see a quantitative description in the annex 1).

Perform linear interpolation of the gray values in the following way: one of the base slice profiles is shifted (see blue profile in [Figure 3], shifted profile in blue dashed line), so that the reference points are congruent; then the congruent profile values of the two profiles along the radial beam are linearly interpolated taking interpolation factor and targeted position into account (result is dashed orange line in [Figure 3]; see a quantitative description in the annex 1).

Position the new radial profile along the radial beam corresponding to the new bone onset point (see final position as orange curve in Figure 3).

Overlay it to the new slice and modify the corresponding gray-values of the voxels touched along this radial ray.

The result is a new cross-sectional structure inside the selected circle.

Typical results for a given radial ray in a linear interpolation are shown in [Figure 4]. One can clearly see the linear morphing of the profile in position and shape into the three new profiles.

Figure 4

typical radial profiles of a linear interpolation approach with two base slices (green and blue, distance 10 mm) and resulting interpolated profiles (IP1 to IP3) of a human tibia shaft near the proximal epiphysis performed with an interpolation factor of 4, therefore three interpolated profiles; threshold value for bone onset detection was set to 100.

In order to avoid sudden transitions between the interpolation sections it is useful to use spline interpolation instead of linear interpolation. We have done this with a cubic spline approach. The processing for the cubic spline interpolation is analogous to the linear spline except that more than two base slices are required, in our case four. Here, the procedure is as following:

Select four slices of equal distance as base slices for the interpolation. The interpolation zone is between the innermost two slices.

Provide the original distance between the base slices and the desired distance between the interpolated slices. the interpolation factor will be calculated as for the linear interpolation case.

Manually: Select a center point in one of the inner base slices.

Manually: Select a maximum radius length (selecting a point outside the bone cross section). Keep care that the corresponding circle is large enough so that it includes all four base cross-sections.

Take one of the inner base slices as the baseline of the new interpolated slice. Do a loop over 360 degrees with an arbitrary angle step size (normally 1° or 0.5°) around the center point and for every angle (and corresponding radial ray) do the following steps:

Generate the four corresponding radial profiles for the base slices.

For every profile, define the onset point (index of profile vector) of bone related gray-value (the outer edge of the bone); these points will be reference points for the following interpolation processes.

Cubic spline interpolation of bone onset: the bone onset for the interpolated profile is interpolated as a cubic spline between the reference points of the two inner profiles for the new target slice taking interpolation factor and targeted position into account (see a quantitative description in the annex 1).

Cubic spline interpolation of the gray values: shift all profiles with respect to one selected so that the reference points are congruent; then the congruent profile values of the two inner profiles along the radius are interpolated as cubic spline based on the values of all profiles and taking interpolation factor and targeted position into account (see a quantitative description in the annex 1).

Position the new radius profile along the radius beam corresponding to the new bone onset point.

Overlay it to the new slice and modify the corresponding gray-values of the voxels touched along this radius.

Our software tool performs most of these actions automatically. The user has to provide a set of coarsely scanned equidistant slices, the interpolation mode (linear or cubic spline) and to specify the desired distance for the interpolated slices. The tool interpolates between the two first scans (the first gap) linearly in all cases because the required lower scan is missing. When cubic spline mode is selected, it employs this mode as long as the required four adjacent base slices are available. For every interpolation gap the tool asks for the selection of a center point and for the selection of a maximum radius point (the issues that are indicated as manualy in the procedure). That can be done graphically by two mouse clicks setting two crosses like in in [Figure 3]. The final gap can, same as the first, only be interpolated linearly. Note that the method is applicable for geometries with approximately circular cross-sections only. The method works best when a closed cortical shell is available.

Typical results for a given radial cubic spline based interpolation procedure are given in [Figure 5]. Depending on the properties of the outer baseline slices the morphing of the new family of profiles (IP1 to IP3) is no more strictly linear in position and shape. For nonlinear geometrical structures for instance near the epiphyses, interpolation with the cubic spline approach yielded generally better results. Therefore, the remainder of this manuscript concentrates on results obtained using the cubic spline method. However, the interpolation between the outermost slices can only be linear because the additional slice for the cubic spline calculation is missing. The difference in computational power is significant (2 seconds for a single linear interpolation process vs. 60 seconds for the cubic spline interpolation on a normal desktop computer using MATLAB R2012b, The MathWorks Inc.), but otherwise there is no real argument to prefer the linear interpolation.

Figure 5

typical radial profiles of a cubic interpolation approach with an interpolation factor of 4; for cubic spline interpolation four base slices (red, green, blue and cyan; distance 10 mm); interpolation factor of 4 provides three interpolated slices (IP1 to IP3) between the inner initial slices (initial lowa, green and initial higha, blue); threshold value for bone onset detection was set to 100; for the cubic spline algorithm the two outermost profiles lowb and highb are required; the boundary conditions are the same as for Figure 4, so that a direct comparison can be made with the linear interpolation result.

A critical element in both the linear and the cubic spline interpolation is the proper detection of the bone onset points. In the algorithm a threshold value can be specified. The spongiosa parts of the epiphyses in particular often have a very thin outer cortical layer that tends to disappear in the pQCT image because of the partial volume effect. The bone seems to be “open” and we only find a weak transition between non-bone gray-value of the outside and the low mean gray-value of the spongiosa. When the algorithm does not identify this point properly it produces unreliable results for the corresponding radius. This situation can normally be detected by the algorithm and warnings can be issued to the operator, so that the threshold value can be optimized.

Quite frequently, the new profile does not cover the entire range between the (new) reference points because of the incomplete overlay of the different, relatively shifted profiles. That means that the innermost voxels around the center point remain unmodified. However, this is not critical when interpolating long bone scans.

Evaluation results

For evaluation of the new method, the proximal 24 cm of a human tibia bone (ex vivo) have been scanned by pQCT from midshaft (position 0) to the proximal tibia plateau (position 240), using a Stratec XCT3000 device (Stratec Medizintechnik, Pforzheim, Germany). The inter-slice distance was 1 mm, and the slice thickness provided by the machine is 2.2 mm. The in-slice resolution was 0.5 mm. The total dimension of a single scan was 655x655 voxels. For the interpolation slices of 10 mm distance have been selected as the baseline cross-sections. Then, cubic spline interpolation with an interpolation factor of 10 has been performed resulting in the reproduction of the nine missing slices at 1mm inter-slice distance. The interpolated cross-sections could now directly be compared to the measured slices.

[Figure 6] shows interpolation results for the region between positions 121 and 131, which are located in the neat shaft area. In this area the changes in shape and cross sectional area are relatively small, so that the interpolation algorithm is not really challenged. By contrast, in the epiphyseal region, which is shown in [Figure 7] (between 194 and 204 mm), there is substantial variation in cross-sectional shape, and the corticalis tends to fade away. As already mentioned this makes interpolation more challenging for any algorithm in epi- and metaphyseal regions than in the shaft. Even so, the optical inspection shows a good adaptation of the outer shape. While most of the contours of the interpolated slices are smooth, we have some minor cripple at the outer edge at the upper right. That is due to the somewhat low quality of the baseline slice 204 at this edge. For our purpose of generating 3D geometric models for FE analysis this is not a critical flaw.

Figure 6

Typical sequence of cross section gray values in the diaphyseal area: top line measured cross sections at positions 121, 124, 128, 131 (in mm); bottom line corresponding interpolated (spline) cross sections 124 and 128, interpolated from base slices 121 and 131.

Figure 7

Typical sequence of measured and interpolated slices in gray values (taken from a representation by the software AMIRA (AMIRA v5.3.3. by Konrad-Zuse-Informationszentrum Berlin + Visage Imaging); upper row measured slices at positions 194, 197,201 and 204 (in mm); the slices at position 194 and 204 have been taken as base slices for a cubic spline interpolation; the resulting slices at positions 197 and 201 are given in the bottom row.

For a better localization [Figure 8] shows an axial cross section of this region.

Figure 8

Axial cross section of the metaphysis area: picture A shows a cross-sectional view when the 10 mm distant base slices (positions 164,174,184,194,204,214) are given only; picture B shows the same view when all measured slices are taken into account (1 mm distance); picture C shows the result, when all linearly interpolated cross-sections are taken; picture D shows the case based on spline interpolation; all views generated by AMIRA.

In [Figure 8] the strengths and weaknesses of our method become quite obvious. The 10 mm resolution leads to the terraced axial cross-section and causes big problems for tools like AMIRA, MIMICS or AVISO when generating 3D geometries out of those slices. The interpolated set of slices can smooth these edges very efficiently. Yet, we can also see that the method did not capture the convex outer envelope at the right upper side. The differences between linear and spline interpolation are not very striking, the left edge being a bit smoother but not a significantly better reproduction of the template.

Moreover, the strengths and weaknesses of our interpolating algorithm become obvious in [Figure 9]. An acceptable smoothing result is obtained in those areas where the bone surface is closed. In the open and topological less distinct area, however, and thus for most of the proximal epiphysis the algorithm fails. It can also clearly be seen, that AMIRA especially has problems in smoothing the distal epiphyseal area where our algorithm shows very good results.

Figure 9

Volume rendering results by AMIRA from left: from full data set (A), from coarse reduced data set (B), from interpolated data set (C); the right picture (D) shows a photo of the scanned bone from a similar perspective.

In order to be able to compare interpolated and measured slices quantitatively we introduced a diversity index as defined in equation 1:

where Hmi is the gray value of voxel i at base slice m and n the number of voxels. As one can see from the equation, our diversity index is the root mean square error (RMSE) divided by a maximally possible RMSEmax, when the main structures of a slice are not overlapping at all. Comparable to the usually given RMSE[7] this index emphasises higher errors. ID is 0 when two slices are identical and it is 1 if the slices are different without any local overlap of the structures inside.

[Figure 10] depicts the variation of the diversity index (ID) with increasing distance from a base slice (position 0). For the epiphyseal region, ID steadily increases and reaches values close to 1 in a distance of 10 mm (dark blue curve). At a distance of only 2 mm, ID reaches a value of 0.45. This reflects the diversity of adjoining pQCT slices in the epiphyseal region. Accordingly, when comparing the reference slice to the interpolated slices in increasing distance (dark red curve) we find nearly the same behaviour, indicating validity of our interpolation algorithm. For the shaft region, diversity from 0 is much smaller with a maximal ID of 0.5, both for the measured as well as for the interpolated slices, which reflects greater preservation of bone ‘shape’ per unit length in the shaft area.

Figure 10

Diversity analysis of measured and interpolated slices of two segments between #194 and #204 (metaphyseal region) and between #111 and #121 (shaft region); presented is either a comparison between a measured or interpolated slice with respect to the lower base slice (given as #0 in the diagram) or a direct comparison between measured and interpolated slice.

When regarding the direct comparison between interpolated and measured slices, given in [Figure 10] as well, we would therefore rank the ID values obtained for diaphyseal interpolation (orange line), which were always well below 0.2, as an excellent result of our interpolation algorithm. For the epiphysis, values around 0.4 suggest that the algorithm is reasonably accurate (green line). A limitation of the diversity index introduced here is the fact that translational shift is not taken into account. That means that completely identical patterns, which simply have been shifted in the slice plane, will be found as different and therefore feature a high diversity index. For that reason, the ID index tends to overestimate, rather than to under-estimate the diversity to a certain extent.

An important quantity for the assessment of biomechanical strength of a bone is the cross-sectional moment of inertia, also known as second moment of area. We calculated a mineral-density-weighted version. Given a bending axis through the center of area of a cross-sectional plane it calculates to

CSMI = Σi d(i)2*A(i)*ρ(i) (2)

where d(i) is the distance of a bone voxel i from the axis, A(i) the cross-sectional area of this voxel and ρ(i) its mineral density.

In [Figure 11] the density-weighted CSMI for two orthogonal bending directions are given along the entire measured length of our tibia model and are compared between interpolated and measured geometries. There is generally good agreement between the interpolated and the measured values, except for the region beyond #204. In that region, we are in the epiphysial regime where the power of our method is limited, as mentioned above. In addition we have differences around #180 and #190 mainly in the y-directed case. Looking for an explanation we found that there were holes drilled in our tibia that reduced the CSMI of the measured slices at these points and in this direction significantly. This of course could not be reproduced by the interpolation. The cavity at #190 seems to influence the base point at #184 which explains the underestimation at #182.

Figure 11

Density-weighted Cross-Sectional Moment of Inertia (CSMI) for the region between #104 to #234. Two different directions of bending have been analysed with x and y indicating their direction of the bending axis in the x-y-plane. Both scanned and interpolated geometries are given. The dotted lines give the position of the base points for interpolation.

In order to demonstrate the suitability of the new interpolation method in combination with a segmentation tool, we used the software package MIMICS (MIMICS Research v 18.0, Materialize) to perform 3D model construction with the different data sets. The result is shown in [Figure 12].

Figure 12

3D-geometries generated by MIMICS from the different data sets: upper row isometric view; lower row coronal cut through the geometry; left side result from the measured high resolution set with 1 mm slice distance; in the middle reduced data set with scans of 10 mm distance; on the right interpolated data set with slice distance 1 mm as well.

The terraced geometry that caused the problems described in the introduction can now be substituted by the 3D-structure that is based on the interpolated data stack and comes quite close to the model generated from the high resolution data set. Yet, there are some differences at the outermost epiphysis. In this region the interpolation mechanism begins to fail because of the lack of clear outer borders of the bone. The surface is a bit more rough, which, however, should not cause problems for use in FE simulations.

Discussion

A powerful method has been introduced that allows for interpolation and generation of intermediate slices between series of pQCT-scans that are too distantly spaced to allow for a generation of a reasonable 3D model. The method has been developed for the 3D-model generation of long bones like the tibia and is therefore adapted for that need. It has been shown that the method reproduces the evolution of the outer and inner shape very accurately, as well as the gray value distribution of the cross sections. This is especially useful for improving the input for segmentation tools like AMIRA, MIMICS or AVISO. The benefit has been shown by comparing results of 3D-model generation with the MIMICS software of the high resolution measured scan set as well as the reduced set and the set containing the high resolution slices recovered by interpolation. It has been found that the quality of the 3D-model based on the interpolated data set is close to a real, highly resolved scan set and significantly better than the results that are achieved when relying on smoothing algorithms in the segmentation tool. Depending on the divergence of shape also larger gaps (for instance in the shaft area) can be reasonably filled with interpolated slices.

As the tool is focused on the processing of long bone scan data it requires a compact structure that shows nearly round cross-sections without substantial elongation or any deep concave bays that go deeper than its center of area. Another requirement is that there should be clear detectable bone outer edges. At the epiphysis our tibia showed spongy-like openings that lead to failures in that area even when using a very low threshold value of 130. So we suggest to limit the application of our algorithm to the diaphysis and metaphysis regions of long bones where we can find a distinct corticalis.

The method can as well be applied for scans from medical imaging sources other than pQCT, e.g. MRI.

We have successfully applied the method in the context of novel exploitation of data from the MUST-study[12], using FE simulation (publication in preparation). In order to perform that analysis, we had to overcome exactly the situation of sparse pQCT scans described in the introduction. Without having a real alternative the interpolated scans allowed us to generate 3D-geometry models of good quality that rendered this analysis possible.

Outlook

At the moment the method is used by us for deriving 3D-models for FE-simulations based on pQCT-scans of studies in the past. The software is currently working on a semi-automatic basis that requires definition of the center point and definition of the working radius. In a next step the implementation of a full automation is foreseen.