R&D LOGO

Journal of Rehabilitation Research and Development
Vol. 39 No. 5, September/October 2002
Pages 609-614

Processing computer tomography bone data for prosthetic finite element modeling: A technical note
Rakesh Saxena, PhD; Santosh G. Zachariah, PhD; Joan E. Sanders, PhD
Department of Bioengineering, University of Washington, Seattle, WA
Abstract—A software scheme is presented to extract the shapes of tibiae and fibulae from amputee computer tomography (CT) data for use in prosthetic finite element modeling. A snake algorithm is implemented to overcome challenges of bone-soft tissue edge detection common in this application. Means to enhance initial guess contours, ensure contour continuity, overcome point-clustering problems, and handle high-curvature regions are also described. Effectiveness of the algorithm is demonstrated on image data from a unilateral transtibial amputee subject.
Key words: active contour model, amputee, fibula, snakes, tibia.
INTRODUCTION

An important challenge in the development of prosthetic finite element (FE) models to study residual limb/prosthetic socket mechanical interactions is accurate specification of the bone shapes within the residu-al limb [1–3]. Several methods can be used to obtain bone shapes, for example, computer tomography (CT),

This material was based on work supported by the National Institutes of Health (HD-31445).
Address all correspondence and requests for reprints to JE Sanders, PhD, Associate Professor; Department of Bioengineering, Harris 309, Box 357962, University of Washington, Seattle, WA 98195; 206-221-5872; fax: 206-221-5874; email: jsanders@u.washington.edu.

magnetic resonance imaging (MRI), and ultrasound. Ultrasonic imaging systems for prosthetics require long scan times and thus are susceptible to subject movement as well as soft tissue distortion because of residual limb submersion in a water tank or mechanical pressure from the transducer [4]. CT does not have these problems and has advantages over MRI in that image quality of the bone, the component of interest here, can be more clearly distinguished.

Despite the capabilities of CT imaging, it is difficult to develop automated or semiautomated methods to extract bone shapes from the CT data for prosthetic FE model applications. Typically, an edge-detection method is used with some form of contouring. However, this approach has two major difficulties. First, the threshold value to achieve effective bone edge detection varies from slice to slice, leading to significant errors [5]. Second, often the cortical bone is so thin, particularly in the proximal tibia, that the edge is not discernable between the soft tissue and the cancellous bone, causing automatic edge-detection algorithms to drift away into the body of the bone. Contouring requires pixel-to-pixel continuity of the edge outline; thus, even sophisticated edge-detection filters have gaps in the outlines that must be rectified manually. Such manual manipulation of the images is tedious and time-consuming.

The purpose of this technical note is to describe an approach for overcoming these two problems. A model-based segmentation method is used, whereby rather than each contour being fixed manually, the model automatically bridges gaps in the contour. So-called “active contour models,” recently applied to nonamputee prosthetic applications in biomedicine, are influenced by global feedback in addition to image features at the local level [6–8]. Thus, they have the potential to overcome the difficulties just described. Means to apply this approach to lower-limb bone extraction are described for use in prosthetic FE modeling application.

THEORY

“Energy-minimizing active contour models” (snakes) were introduced by Kass et al. in 1988 [9]. The term “active” refers to the deformability or flexibility of the snake to minimize its energy dynamically in order to adapt it to given image data. In the present application, snakes are used to identify the boundary between bone and soft tissues in the lower limb. The boundary extraction process is initiated by placing an estimated contour close to the intended boundary within the image. The snake algorithm then adjusts this estimate by minimizing the snake’s total energy so that it locks onto the boundary. The total energy includes a feature energy term and additional terms to ensure contour smoothness. The total energy functional of the snake, represented by a contour v(s) = (x(s),y(s)), where s is the arc length, can be written as [9]

equation 1.

Eint is the internal energy due to forces of bending and stretching that maintain snake smoothness. Eimage is the energy due to image forces that push the snake to take the shape of the image features of interest. In the present application, we used the two-dimensional gradient of pixel intensity to determine image energy, since the boundaries appear as local maxima in the gradient images. Eint can be written as

equation 2

where the subscript s denotes the derivative. In equation (2), the first term accounts for stretching of the snake and will have a large value where there is a gap in the curve. The second-order term accounts for bending and will be large in a region of high curvature. Thus, the values of a and b locally control the stretching and bending of the snake at that point along the contour.

Amini et al. proposed a dynamic programming algorithm to discretize the variational problem given by equation (1) [10]. Using the finite difference method, one can approximate the two derivatives in equation (2) as

ferivatives of equation 2

formula 3

and

result of formula 3

formula 4

The total snake energy is given by

formula for total snake energy

where the summation is done over all the points. The implementation by Bregler and Slaney of this algorithm in MATLAB (MathWorks, Inc., Natick, Massachusetts) provided the starting point for the present study [11].

METHODOLOGY AND RESULTS

Application of the method is illustrated as follows. The residual limb of a unilateral, lower-limb, male amputee aged 55 was imaged with the use of a General Electric CT scanner (Genesis Hispeed RP). The subject was diabetic and underwent a limb amputation 2 years prior because of an industrial injury. The residual limb was positioned within a patellar-tendon-bearing socket during CT scanning. Scanner settings were 120 kV, 120 mA, with a slice thickness of 3 mm and a slice spacing of 3 mm.

Each image slice was processed first with the adjustment of the upper and lower bounds of the grayscale image so that the air and some soft tissue could be masked, while retaining good contrast at the boundary between bone and soft tissue. This procedure helped to highlight the features of interest.

Next, snake parameter settings were set so that the contribution of the three energy terms (stretching, bending, and image) were within the same order of magnitude. The two principal parameters are (1) acontrols the energy contribution by the gap between adjacent boundary points and (2) b—controls the energy contribution caused by the bending of the curve. In a closed contour, increasing a has the effect of pushing points toward equal spacing along the contour. In this application, a is set to zero so that snake points are not penalized for moving toward strong edges or for clustering.

Additional settings affect the resolution of the contour: Dx and Dy - the maximum number of pixels moved in the x and y directions, respectively, which determines the size of the neighborhood for the minimum energy search, and ex and ey - the pixel resolution in the x and y directions, respectively. For the algorithm to run, each snake point is moved in the neighborhood formula for snake poing meighborhood, and the total minimum energy given by equation (5) is calculated at the end of every iteration. The iterations are continued until a global minimum in total energy is reached.

For each slice, the determined contour depends on an appropriate initial guess. To overcome initial guess problems in slices below the most proximal slice, the final snake boundary from the previous slice is used as the initial boundary guess for the next distal slice. Because the bones tend to reduce in size from proximal to distal, the initial guess contour is usually larger than the actual bone boundary. During iterations of the algorithm, the snake is pushed toward the actual boundary.

A potential limitation specific to this prosthetics application is that because equations (1) to (5) apply to open snakes (i.e., the first point and the last point of the snake are not constrained to overlap), unfolding can occur. This unfolding is due to contributions from the a and b terms in equation (2) trying to open up and shorten the boundary to reduce the bending and stretching energies. Figure 1(a) and (b) is a CT slice illustrating this difficulty. Note that in equations (3) and (4), the bending and stretching energies are not completely defined at the ends of the contour. Equations (3) and (4) are therefore modified to ensure continuity of the internal energy of the snake over a closed contour, by defining

internal energy formula

equation 6.

Figure 1.
Contour opening and point clustering. (a) and (b) Contour of tibia opening and (b) a magnification of (a). (c) Closing of contour with use of techniques described to ensure continuity. (d) Result of cubic spline fitting and resampling to eliminate point clustering.

Figure 1(c) shows the new image after implementing this change. When the feature energy is uniform, such a closed snake will converge to a circle with a radius where the bending energy exactly matches the stretching energy.

Another challenge in application to prosthetics is that in local regions of high image gradient values and/or high snake curvature, the snake points tend to move closer together, thereby decreasing the total snake energy. Rather than being a disadvantage, this is an advantage in the present application, since the initial guess moves away from points where the edge is weak as well as provides better contour definition in high curvature areas. However, since the final solution for one image is used as the initial guess for the next image in the series, this clustering tends to accumulate over several slices. For this limitation to be avoided, a cubic spline is fitted to each solution and then the curve is resampled to generate equally spaced points for the next guess. Resampling also permits the number of snake points to be reduced as the total contour length decreases toward the distal images. The result of cubic spline fitting and resampling of Figure 1(c) is shown in Figure 1(d).

A further complication in application to prosthetics is that in regions of high curvature where the bending energy dominates over the feature energy, the snake tends to cut off contour sections of large curvature (Figure 2(a)). To overcome this challenge, we adjusted the local values of b(s) so that local bending energy was constant at all points on the initial guess for an edge contour for a given slice. This ensured that, in the absence of other forces, the snake tends toward the shape of the current guess rather than toward a circle. The local values of b were determined so that the global b (total bending energy divided by total feature energy) retained the desired value—0.26 in the present case. This provided good overall results (Figure 2(b)), ensuring that the bending energy exerts a uniform push-pull on the snake, irrespective of the magnitude of the feature energy and the local curvature.

Figure 2. cutting off of contours.

Figure 2.
Cutting off of contours in regions of high curvature. (a) Use of a constant b over entire contour. (b) Use of a variable b but constant bending energy over entire contour. A variable b better follows the fibular boundary.

Once all images are processed, points on the contours establish the bone surface. Results for the tibia of a trans-tibial amputee are shown in Figure 3.

Figure 3. Reconstruction of bone shape.

Figure 3.
Reconstruction of bone shape. Results from tibia of a transtibial amputee subject are shown. Units are in millimeters.

DISCUSSION

The methods that we presented using active contour modeling help to overcome challenges in extracting tibia and fibula bone shapes from CT image data of transtibial amputee subjects. Difficulties with threshold level varia-bility from slice to slice are overcome, as are challenges in defining the bone and/or soft-tissue boundary in regions where the cortical bone is very thin or absent. Means for establishing initial guess contours, ensuring continuity, overcoming point clustering problems, and handling high-curvature regions in the images facilitate effective performance of the algorithm for use specifically for prosthetic FE modeling efforts.

REFERENCES
1. Silver-Thorn MB, Steege JW, Childress DS. A review of prosthetic interface stress investigations. J Rehabil Res Dev 1996;33:253–66.
2. Zachariah SG, Sanders JE. Interface mechanics in lower-limb external prosthetics: A review of finite element models. IEEE Trans Rehabil Eng 1996;4:288–302.
3. Zhang M, Mak AF, Roberts VC. Finite element modelling of a residual lower-limb in a prosthetic socket: A survey of the development in the first decade. Med Eng Phys 1998;20:360–73.
4. He P, Xue K, Chen Q, Murka P, Schall S. A PC-based ultrasonic data acquisition system for computer-aided prosthetic socket design. IEEE Trans Rehabil Eng 1996;4:114–19.
5. Sumner DR, Olson CL, Freeman PM, Lobick JJ, Andriacchi TP. Computed tomographic measurement of cortical bone geometry. J Biomech 1989;22:649–53.
6. Sageflat JH, Olstad B, Torp AH. Automatic custom femoral stem design based on active contours. Proc SPIE, Med Imaging Image Processing 1995;2434:426–37.
7. Snel JG, Venema HW, Grimbergen CA. Detection of the carpal bone contours from 3-D MR images of the wrist using a planar radial scale-space snake. IEEE Trans Med Imaging 1998;17:1063–72.
8. Safont LV, Marroquin EM. 3D reconstruction of third proximal femur (31-) with active contours. Proceedings 10th International Conference on Image Analysis and Processing; 1999; Venice, Italy. Los Alamitos (CA): IEEE Computer Society; 2002. p. 458–63.
9. Kass M, Witkin A, Terzopoulos D. Snakes: Active contour models. Int J Comput Vis 1988;1:321–31.
10. Amini AA, Weymouth TE, Jain RC. Using dynamic programming for solving variational problems in vision. IEEE Trans Pattern Anal Machine Intell 1990;12:855–67.
11. Bregler C, Slaney M. SNAKES, Interval Research Corporation; 1995.
Submitted for publication August 3, 2001. Accepted in revised form May 2, 2002.