f

Modelling and simulation of the intervertebral movements of the lumbar spine using an inverse kinematic algorithm L. W . S u n I

R.Y.W.

Lee I

W . Lu z

K.D.K.

Luk z

1Department of Rehabilitation Sciences, The Hang Kang Polytechnic University, Hang Kang 2Department of Orthopaedic Surgery, The University of Hang Kang, Hang Kang

A b s t r a c t - - A n inverse kinematic m o d e / i s presented that was employed to determine the o p t i m u m intervertebral j o i n t configuration for a given forward-bending posture of the human trunk. The lumbar spine was modelled as an open-end, kinematic chain of five links that represented the five vertebrae (L 1-L5). An optimisation equation with physiological constraints was employed to determine the intervertebral j o i n t configuration. Intervertebral movements were measured from sagittal X-ray films of 22 subjects. The mean difference between the X-ray measurements of intervertebral rotations in the sagittal plane and the values predicted by the kinematic model was less than 1.6 °. Pearson product-moment correlation R was used to measure the relationship between the measured and predicted values. The R-values were found to be high, ranging from 0.83 to 0.97, for prediction of intervertebral rotation, but poor for intervertebral translation (R-0.08-0.67). It is concluded that the inverse kinematic model will be clinically useful for predicting intervertebral rotation when X-ray or invasive measurements are undesirable. It will also be useful to biomechanical modelling, which requires accurate kinematic information as m o d e l input data.

Keywords--lnverse kinematics, Optimisation, Intervertebral motion, Lumbar spine

\

Med. Biol. Eng. Comput., 2004, 42, 740-746

1 Introduction

INSTABILITYOF the human spine is a significant cause of low back pain, which has resulted in substantial financial costs to the health services, in Britain, this amounts to approximately £650million per annum (ROSEN et al., 1994). Knowledge of the intervertebral movements of the lumbar spine is clinically useful in the assessment of spinal disorders (such as instability) and treatment outcome. Measurement of the intervertebral movements of the lumbar spine is fraught with difficulties. This is primarily because the spine is rather inaccessible, and the nature of the movements is very complex. Radiographic, electro-optical and electromagnetic techniques have been used in the measurement of lumbar spine motion with some success (BROWN et al., 1976; PEARCY, 1985; PEARCY and HINDLE, 1989; DVORAK et al., 1991; LEE, 1995; 2001; LEE and EVANS, 1997; LEE and MUNN, 2000; TAKAYANAGI et al., 2001). However, radiographic techniques are complicated and have the inherent health risk of repeated X-ray exposure.

Correspondence should be addressed to Dr Raymond Y. W. Lee; emaih [email protected] Paper received 5 February 2004 and in final form 1 June 2004 MBEC online number: 20043934 © IFMBE: 2004 740

J

Surface measurements using markers or sensing devices are subject to large error owing to the deformation of underlying soft tissues disguising the true vertebral movement (PEARCYand HINDLE, 1989). Such techniques are only accurate and reliable for determining the total movement in regions of the spine, as opposed to discriminating the contribution of the movement of individual intervertebral j oints (PEARCYand HINDLE,1989; LEE, 2001). Surface measurements of intervertebral movements can only be achieved by fixing markers or sensing devices rigidly to the spinous processes with pins or screws (KAIGLEet al., 1992; KANAYAMA et al., 1996), but this is an invasive surgical procedure. The human spine has redundant degrees of freedom, giving us great flexibility in the performance of a movement. For instance, a forward bending movement of the spine can be accomplished with an infinite number of combinations of configurations of the various intervertebral joints, inverse kinematics is a method used in robotic engineering for determining joint configurations given a desired position and orientation of the end effector of the robot in achieving a certain goal (CRAIG, 1989; MCCARTHY, 1990; ALLARD et al., 1995; ZHANGand CHAFFIN,1996; ZATSIORSKY, 1998). We propose that an inverse kinematic algorithm can be used to determine whether a certain combination of intervertebral joint configurations may be the most feasible and reasonable way of performing forward bending. The inverse kinematic method has been successfully used to model the movements of various body joints (SOMMERand Medical & Biological Engineering & Computing 2004, Vol. 42

MILLER,1980; FUJIE et aL, 1993; JUNGet aL, 1994; ZHANGand CHAFFIN,1996; ZATSIORSKY,1998; WANG 1999). For instance, JUNG et al. (1994) employed a psychophysical cost function (potential function) to define a cost value for each joint movement angle, and developed a regression model to predict the perceived discomfort with respect to the joint movement. ZHANG and CHAEEIN (1996) developed a new optimisationbased differential inverse kinematics approach for modelling three-dimensional dynamic seated postures. WANG (1999) also utilised an inverse kinematics algorithm to predict arm prehension postures. However, the inverse kinematic method has never been applied in the case of the spine. This is probably because the spine has an unlimited number of degrees of freedom, and finding an optimum solution to the inverse kinematic problem would be extremely difficult. However, a solution could be feasible if appropriate constraints could be imposed on the kinematic modelling. For instance, the positions of the most posterior parts of the spinous processes can be determined by surface measurements. Such information can be used to reduce the number of solutions. In addition, an optimisation function can be employed to determine the best solution. The purpose of this study was to evaluate the validity of using an inverse kinematic model to determine the intervertebral joint movements, with given knowledge of the total flexion movement of the lumbar spine and the positions of the spinous processes.

L1

6~

L2

e~

d2C ~ ~

xs

L3 s2

L4

L5 s

sacrum

~

x

2 Methodology 2.1 Prediction o f intervertebral movements using an inverse Idnematic algorithm

The lumbar spine was modelled as a five-link system from the L1 to L5 vertebra (Fig. 1). it was assumed that there were 3 degrees of freedom (DOFs) (one rotation and two translations) for each intervertebral joint. Forward flexion of the spine was assumed to be confined to the sagittal plane (PEARCY, 1985). The Denavit-Hartenberg convention (McCARTHY, 1990; ZATSIORSKY, 1998; TAKAYANAGI et al., 2001) was used to describe the multilink chain. The anterosuperior comer of the sacrum was considered to be the origin of the kinematic chain in the global co-ordinate system. A local co-ordinate system was defined for each vertebra, with the x-axis running from the anterosuperior comer to the posterosuperior comer of the vertebrae, and the y-axis being perpendicular to the x-axis (TAKAYANAGIet al., 2001). The following information was assumed to be known: (a) The positions of the most posterior parts of the spinous processes: in this study, the positions were obtained from radiographs. However, clinically, these bony landmarks can be easily palpated through the skin surface. Their positions can be predicted from surface measurements (STOKES et al., 1987) or directly measured by noninvasive techniques, such as ultrasonic scanning (PORTER et al., 1978; SUZUKI et al., 1989). (b) The total movement of the whole lumbar spine: this was determined from the changes in the curvature (lordosis) of the lumbar spine during the forward bending movement. Lumbar lordosis was defined by the angle of intersection between a line nmning along the inferior border of T12 and a line along the superior border of the sacrum. Clinically, this was referred as to as the Cobb's method and could be readily determined by surface measurement techniques, (c) The geometry of the vertebra and the length of each link of the kinematic chain: these parameters were obtained Medical & Biological Engineering & Computing 2004, Vol. 42

Fig. 1 Five-line kinematic model of lumbar spine. Four corners" o f each vertebra (A,B,C,D) are identified by fitting vertebral body image with quadrangle. Local co-ordinate system o f intervertebral joint (shown for L2/3 segment only) is" defined by xcaxis joining posterosuperior and anterosuperior corners" o f lower vertebra o f joint, with origin at point D. yi-axis was perpendicular to xi-axis'. Position of spinous process is" denoted as" Si. di represents length o f kinematic chain from ith vertebra to (i + 1)th vertebra. Global co-ordinate is"shown as x- and y-axes

from previous studies (TENCER and MAYER, 1983; PANJABI et al., 1992b; CHAFFIN et al., 1999). The inverse kinematic problem was to derive the intervertebral joint configuration with knowledge of the positions of the spinous processes. The joint angles Oi and the x andy translations of the vertebrae (LI-L5) were the state variables COl.Thus there were three unknown variables for each joint, and, with five joints in the kinematic chain, the total number of unknown variables was 15. The co-ordinates of the most posterior parts of the spinous processes Si were expressed as functions of the state variable coi, which was a function of joint angles Oi and the xi and Yi translations of the ith vertebrae (i = 1, 2 . . . 5) Si = f ( ~ i )

(1)

where ~ i = [Oi, xi,Yi] ~

741

To derive the variable coi with a given S/, the inverse of the equation would be required, it denotes

coi = f

l(Si)

(2)

Solving the above equation would be difficult. There was more than one solution for a given set of S values. The inverse kinematic problem was solved by the general equation & -- J + - S + (I - J + - J )

-k G

(3)

where J was the Jacobian matrix J(~) = ~

3

(4)

s

J + was the psuedo-inverse of J, and I was the identity matrix. The term P, referred to as the potential function, defined secondary conditions that needed to be fulfilled, it eliminated redundancy, accounted for the physiological limits of the intervertebral joints and avoided kinematic singularities. An optimum solution of the inverse kinematic problem was obtained by minimising the potential function (VAN DEN BOGERT et al., 1994; TILG et al., 2000). The following potential ftmctions were employed in this study: (i)

The error e in predicting the total movement of the lumbar spine 70 was minimised. The potential ftmction was /'1(~)

= 1 ~2

(5)

where e 70 - 81 - 02 - 83 - 84 - 85. The following equation was then used to obtain values of 8 i for the error to be less than a pre-defined value (1°).

OP1 &o

=(70-01-02-03-04-85)[1

function 1°3, so that an approximate solution of co could be determined. P3(co) = - ~ / d e t ( J - J+) = - w , P --+ Pmi,,D C 0

(9)

An iterative refinement algorithm was used to improve the solution (CHEN et aL, 1999). The positions of the spinous processes were predicted using the approximate solution of co and (1). The residuals r of the solution were defined as the difference between the predicted and actual spinous processes' positions. The solution was accepted if the residual was less than its 2-norm II r II 2. if not, the inverse kinematic problem was considered to be ill-conditioned. Singular value decomposition was employed (GoLus and VAN LOAN, 1996; CHEN et al., 1999) to estimate the error. The approximate solution of co was adjusted accordingly. The new Oi, xi andyi values would then serve as initial values for obtaining another solution of co using the inverse kinematic equation and the potential functions. This process was repeated until the residual of the solution was less than the 2-norm. The inverse kinematic problem was then considered to be well defined, and the solution of co was considered optimum. The inverse kinematic algorithm for obtaining an optimum solution of co is summarised in Fig. 2. 2.2 Experimental validation Lateral radiographs of the lumbosacral spines of 22 subjects (nine men and 13 women, mean age = 40-4-14 years) were obtained from the Duchess of Kent Children's Hospital, HK. They were taken with the lumbar spine in full flexion and extension while the subjects were standing. The radiographs were either taken during routine clinical examination or used in other research studies. Subjects were diagnosed with non-specific

1 1 1 10...0] set initial value of state variable m

(6) (ii) The potential function for constraining the intervertebral rotations and translations (Oi, xi and Yi) within the physiological limits of the joints was 5 i=1 5 p2(co) =

begin iteration loop

J

½ (8i mox - 8i) 2 + (8~ m~, - 8~)

C ½(Xim,x __ Xi) 2 + (Ximi, __ Xi) 2

i=1 5

find Jacobian

4,

(7)

define potential functions and apply constraints

+

Z ½(y, ~°x - y,) ~ + (y, ~,, - yi)

i=1

where 8 i rain, 8i . . . . . x i rain, Xi . . . . . Y i rain, Yi . . . . were the minimum and maximum values of 8i, xi and yi. The minimum and maximum allowable joint angles were 5 ° and 22 ° , respectively. These values were based on experimental results reported by previous authors (PEARCY, 1985; SCHAFFER et al., 1990; PANJABI et al., 1992b). The minimum and maximum allowable values of the x and y-translations were 1 0 m m and 10mm. Previous studies showed that the intervertebral joints did not exhibit translation beyond these values (WHITE and PANJABI, 1978; PANJABI et al., 1992b). it follows that -

+ update values of m ]

obtain approximate solution for optimum intervertebral joint configuration oY

determine residual r

I

(x~ m°x - - x~) 3 + (x~ mi, -- X~) 3

(8)

N ~

y

is r< p-norm?

( y , m.x - - Y*) 3 + (y, m*. - - Y*) 3

(iii) There were a number of possible solutions for 8i, xi and yi after the potential functions P1 and P2 were employed. The Jacobian determinant P3 was then used to eliminate those solutions with kinematic singularities. Minimising the determinant would allow convergence of the potential 742

determine 6)

{ (Simax-Si) 3-]-(Simin-Si) 3 }

3P5 a(o

determine unique pseudo-inverse of Jacobian providing minimum norm solution

final solution for intervertebral joint configuration m

Fig. 2 Inverse kinematic algorithm for predicting intervertebral

movements" of lumbar spine Medical & Biological Engineering & Computing 2004, Vol. 42

Table 1

Mean measurement error (-- 1 SD) and mean intra-class correlation coefficient (ICC) values for radiographic measurements

Mean measurement error

Angle of rotation, deg Postero-anterior x translation, mm Supero-inferior y translation, mm

Table 2

Mean ICC (3,1)

L5/S

L4/5

L3/4

L2/3

L1/2

L5/S

L4/5

L3/4

L2/3

L1/2

0.3 0.02 0.02

0.4 0.03 0.02

0.3 0.03 0.02

0.5 0.06 0.02

0.3 0.03 0.03

0.96 0.99 0.99

0.96 0.99 0.99

0.96 0.99 0.99

0.88 0.98 0.98

0.92 0.95 0.95

Comparison of mean measured and predicted kinematic parameters" of 5 intervertebral joints

Motion segment Angle of rotation, deg Mean measured results from X-ray film i Mean predicted results j Mean absolute difference between measured and predicted results i j Correlation between measured and predicted values

L5/S 8.0 (6.0) 6.9 (6.3) 1.6 (1.2)

L4/5 10.6 (4.8) 10.5 (5.0) 1.1 (0.8)

L3/4

L2/3

L1/2

9.0 (3.5) 9.1 (4.0) 1.0 (0.7)

7.1 (3.3) 6.7 (3.0) 1.0 (0.7)

5.2 (2.6) 4.7 (3.0) 1.2 (0.9)

0.95

0.97

0.87

0.91

0.83

0.17 (0.19) 0.12 (0.16) 0.14 (0.10)

0.31 (0.19) 0.27 (0.23) 0.12 (0.08)

0.27 (0.15) 0.21 (0.19) 0.16 (0.11)

0.22 (0.17) 0.12 (0.21) 0.11 (0.08)

0.20 (0.25) 0.16 (0.31) 0.17 (0.12)

0.15

0.67

0.15

0.34

0.64

0.03 (0.08) 0.02 (0.11) 0.05 (0.04)

0.02 (0.08) 0.04 (0.14) 0.08 (0.05)

0.03 (0.06) 0.03 (0.12) 0.04 (0.03)

0.05 (0.07) 0.05 (0.12) 0.08 (0.05)

0.03 (0.08) 0.01 (0.13) 0.05 (0.04)

0.59

0.08

0.31

0.22

0.64

x-translation, cm

Mean measured results from X-ray film i Mean predicted results j Mean absolute difference between measured and predicted results i j Correlation between measured and predicted values y-translation, cm

Mean measured results from X-ray film i Mean predicted results j Mean absolute difference between measured and predicted results i j Correlation between measured and predicted values

LBP with no pathologies. Subjects who showed any signs of fracture or dislocation, spinal instability, spondylolisthesis, narrowed disc spaces, osteophytes, transitional lumbosacral vertebrae or any structural disorders of the lumbar spine, or those who had previous history of spinal surgery, were excluded. The positions of the vertebrae were identified on the radiograph by the fitting of quadrangles around the vertebral bodies (Fig. 1). The positions of the most posterior parts of the spinous processes were also recorded. Images of the inferior vertebra of a motion segment of the flexion and extension radiographs were superimposed, and the images of the superior vertebra of the segment on the two films were then compared. Intervertebral rotation was determined by the change in the angle of rotation of the superior vertebra. Intervertebral translations along the x- and y-directions were given by the changes in the locations of the centre of the superior vertebra. This measurement method has been described in detail in previous studies (LEE, 1995; LEE and EVANS, 1997). The radiographic measurements of intervertebral movements were repeated five times. Intra-class correlation (ICC) (3, 1) was employed to examine the repeatability of measuring the three kinematic parameters among the five measurements (Table 1). The ICC value ranged from 0.88 to 0.97 for intervertebral rotation, and from 0.95 to 0.99 and from 0.95 to 0.99 for xandy-translations, respectively. The radiographic measurements were considered to be sufficiently reliable and accurate. They were then compared with the intervertebral movements predicted by the inverse kinematic algorithm so that the validity of the algorithm could be established. 2.3 Statistical analysis The Statistical Package for Social Sciences (version 11 for Windows)* was employed in this study. Linear regression was *SPSS Inc., Chicago, Illinois 60606, USA

Medical & Biological Engineering & Computing 2004, Vol. 42

used to examine the relationship between the radiographic measurements and the values predicted by the inverse kinematic algorithm. The Pearson product correlation coefficient was determined to assess the degree of association. This analysis was performed for each of the three kinematic variables, namely, the intervertebral rotation and the intervertebral translations along the x- and y-directions. 2.4 Sensitivity analysis Previous studies showed that there could be significant intersubject variations in vertebral geometry, particularly in spines with degenerative changes (ZHOU et al., 2000). Inaccuracies in the geometry data of the inverse kinematic model could lead to errors in the predicted intervertebral rotations and translations. In addition, the physiological constraints, as defined by (7), were based on normal subjects. The minimum and maximum allowable joint movements might need to be altered in subjects with instability or other pathological diseases. Sensitivity analysis of the model was thus performed by changing the input data by 10% and evaluating the corresponding changes in the predicted intervertebral movements. The following input data were examined: the length of the kinematic chain, the position of the spinous process, and the physiological constraints of (7).

3 Results Table 2 provides the description statistics of the values obtained by radiographic measurement and those predicted by the inverse kinematic algorithm. The prediction error, which was defined by the differences between the measured and predicted values, and the correlation between the two values are also presented in the Table. The values are presented separately for each kinematic parameter and for each intervertebral joint. 743

2t 20 "o

15'-I 10 4

•

,'.~ ."

:, i } ~ " " R=o94

g'~-1'5

-1'0

-5

"~)

g

1'0

1'5

2'0

2'5

,,5 .-_o

measured angle of rotation, deg Fig. 3 Regression between the measured and predicted intervertebral rotations for different motion segments'. (e) L1/L2; ( x ) L2 /L3; (D) L3/L4; (A) L4/L5; (*) L5/5; (-) regression line fit

3.1 Intervertebral rotation

The mean predicted values of intervertebral rotations ranged from 4.7 ° to 10.5 ° (Table 2). The mean predicted error was found to be between 1.0 ° and 1.6 °. Fig. 3 shows the regression between the predicted and measured values. The correlation between the two sets of values was high, with R-values ranging from 0.83 to 0.97 (Table 2). This suggests that the inverse kinematic algorithm can be used reliably to predict the intervertebral rotation. The accuracy of prediction varied with different intervertebral segments, with the L1/2 segment being the least reliable level. 3.2 Intervertebral translation The mean predicted values of the x- and y-translations ranged from 0.12 to 0.27mm and from 0.01 to 0.05 mm, respectively (Table 2). The degree of correlation between the predicted and measured values of the translations was poor, with R ranging from 0.15 to 0.67 for the x-translation and from 0.08 to 0.64 for the y-translation. 3.3 Sensitivity analysis Table 3 shows the mean percentage changes in the predicted intervertebral movements as a result of changes in the various input data. The changes in the translation movements were the largest. The mean changes in the angle of rotation due to 10% changes in kinematic chain length and spinous process position were 10.4 -4- 3.5% and 12.0 -4- 7.1%, respectively. However, the angle of rotation was relatively unaffected by changes in the physiological constraints of (7).

4 Discussion 4.1 Prediction o f intervertebral movements using an inverse Idnematic algorithm

The present study demonstrates that the inverse kinematic algorithm can be used reliably to predict the intervertebral rotations of the lumbar spine segments but not intervertebral

translations. The magnitude of the lumbar intervertebral rotations in the sagittal plane was reported to be 13-16 °, and the magnitude of the intevertebral translations was between 1 and 3 mm (WHITEand PANJABI, 1978; PEARCY and HINDLE, 1989; YAMAMOTO et al., 1989; PANJABI et al., 1992b). The measured and predicted values of intervertebral rotations were in close agreement with those reported in previous work, but the magnitude of intervertebral translation observed in our study was smaller than the values previously reported. The small magnitude of the translational movement might explain why it was extremely difficult or ever impossible, to predict the movement reliably, in addition, the mean total magnitude of sagittal rotation of the whole lumbar spine (from L1/2 to L5/S1) was found to be 39.7 ° in this study. This was similar to that reported in previous studies (WHITE and PANJABI, 1978; PEARCY and HINDLE, 1989). The differences between the measured and predicted values of intervertebral rotations were generally small (less than 1.6 °). The prediction error was found to be high for the L 1/2 segment. This could be because the segment was at the end of the kinematic chain, and error accumulated along the chain. Fig. 3 shows excellent correlation between the predicted and measured values of intervertebral rotations. The established regression equation can be employed further to enhance the accuracy of the prediction. However, it should be pointed out the present equation may only be applicable to middle-aged subjects with non-specific low back pain. Different regression equations may have to be established for other population groups. it is suggested that the inverse kinematic technique can be used to predict intervertebral rotation in a clinical situation when it is not desirable directly to measure the movements using radiographs or other methods with health risks. Another potential application of the results of this study is biomechanical modelling. Past mathematical models of the lumbar spine (MCGILL and NORMAN, 1986) made various assumptions about the kinematics of the spine. For instance, it was suggested that, during flexion, the relative magnitudes of movements in the various motion segments were 13.2% for L1/2, 21% for L2/3, 29% for L3/4 and 23.6% for L4/5. it was also assumed that all lumbar segments moved together simultaneously, and that there was a constant relationship between the angle of rotation at each lumbar segment. in this study, it was shown that the distribution of movements was 13%, 18%, 23% and 27% among the L1/2-L4/5 segments, which generally agrees with the figures provided by MCGILL and NORMAN (1986). However, in furore modelling work, the actual intervertebral rotations could be predicted by the present inverse kinematic model, and no assumptions have to be made regarding the distribution of movements among the segments. The relationships between the various segments did not have to be assumed to be constant throughout the movement. The accuracy of biomechanical models can be tremendously improved. Knowledge of the intervertebral translations of the spine may be useful in the clinical assessment of some spinal pathologies such as spondylolysis or spondylolisthesis. Unfortunately, it appears that this can only be reliably determined by radiographic measurement, and the inverse kinematic algorithm was found to be rather unreliable in predicting such movement. This may be

Table 3 Mean pereentage changes in predicted intervertebral movements" due to 10% changes in input data of inverse kinematic model

Mean percentage changes in predicted intervertebral movements Input data (10% change)

angle of rotation

x-translation

y-translation

Length of kinematic chain Spinous process position Physiological constraints

10.4% (3.5%) 12.0% (7.1%) 7.8% (3.0%)

27.1% (15.0%) 30.4% (23.0%) 22.7% (13.2%)

19.62% (10.8%) 25.3% (16.7%) 18.9% (6.3%)

744

Medical & Biological Engineering & Computing 2004, Vol. 42

an area for furore research. The kinematic model could be further refined to increase the prediction of such movement, and noninvasive methods could be developed to measure the movement directly. The inverse kinematic model developed in this study can be modified to determine the intervertebral movements of other regions of the spine. For instance, the lordosis and spinous process positions of the cervical spine can be measured clinically, so that the kinematic mechanisms of healthy and painful necks can be studied. However, an inverse kinematic model of the cervical spine is likely to be more complex than the present model, as it involves seven motion segments and a larger number of degrees of freedom of movement. 4.2 Potential functions To determine the inverse kinematics of a redundant kinematic chain, a large number of solutions may be possible, and it is necessary to determine the optimum solution using some potential functions. In this study, one of the potential functions employed was the total movement of the lumbar spine. This is an appropriate choice, because the movements of the individual joints should be dependent on the total movements of the spine, it can be argued that to minimise energy expenditure and enhance efficiency of movement, a person should perform forward bending of the trunk with as little movement of each intervertebral joint as possible, and the magnitude of the intervertebral movement should not exceed the physiological limits of the joints. Constraints were thus imposed on the kinematic links, and potential functions were employed to minimise the magnitude of the intervertebral movements. The constraints imposed were based on the physiological range of movements reported in previous studies (WHITEand PANJABI, 1978; PEARCY and HINDLE, 1989; YAMAMOTO et al., 1989; PANJABI et al., 1992b). The Jacobian determinant was also employed as a potential function in deriving the inverse kinematic solution. The design of a kinematic mechanism loses at least one degree of freedom when singularity occurs. This happens when the Jacobian determinant becomes zero. For a redundant kinematic mechanism, its manipulability would be minimum when the Jacobian determinant was minimised. Convergence of this potential function would allow a solution to be determined. The fidelity of the solution was examined after the potential functions had been applied. This was achieved by predicting the positions of the spinous processes using the solution and comparing the predicted positions with the known positions. An iterative procedure was employed to minimise the residuals, thus ensuring the optimum solution was obtained. The 2-norm of the residuals was used as the optimisation criterion as it was most computationally efficient when compared with the 1-norm and infinity-norm (GOLUB and VAN LOAN, 1996; CHEN et al., 1999). Previous researchers have employed other techniques to solve the inverse kinematic problem with redundancy. For instance, as opposed to the algebraic method employed in the present study, WANG (1999) employed a geometric method to predict arm reaching posture. The method was based on previously observed behaviours of the arm and was able to solve the kinematic problem in a straightforward way. The advantage of this method was that no matrix inverse calculation was required, avoiding the stability and convergence problems occurring near a singularity of the Jacobian. However, owing to the complexity and unpredictable nature of spinal motions, previous research has not been able to establish any behavioural rules for intervertebral movements, it is not possible to employ a behaviour-based approach for predicting spinal movements, and the pseudo-inverse method has thus been chosen. The potential functions and the iterative procedure were used in this study to solve the convergence problem. Medical & Biological Engineering & Computing 2004, Vol. 42

4.3 Sources o f error it was shown that there would be small changes in the predicted values ofintervertebral rotation if there were inaccuracies in the geometric data of the vertebrae. In our pilot study, we measured the actual geometry of the vertebral images and compared it with the anthropometric data employed in the model. The differences between the two sets of data were found to be less than 10%. The sensitivity analysis showed that such differences would only lead to small errors in the prediction of the intervertebral rotation. However, it was shown that the predicted values of intervertebral translation were very sensitive to changes in input data. This could be a reason why the prediction of translation was inaccurate. The present model may not be appropriate for patients with significant vertebral deformities. The geometry of the spines of these patients will be grossly different from the anthropometric data obtained in the literature. Geometric data will have to be determined using radiographs or other imaging techniques for accurate model prediction. On the other hand, the sensitivity analysis showed that the predicted intervertebral rotation was insensitive to changes in the physiological constraints of the kinematic model. The model may therefore be employed for patients with hypo- or hyper-mobility of the spine. Although radiographic measurements were used as the 'standards' for comparison with the values predicted by the inverse kinematic algorithm, it should be pointed out that they were not free from errors. The accuracy of radiographic measurements may be affected by the clarity of image, the number and positions of chosen landmarks or markers, the process of tracing and superimposition, radiographic quality, intra- and inter-observer variance, the measurement method and the magnitude of the measured motion (PEARCYand HINDLE, 1989; PANJABI,1992a; LEE and EVANS, 1997). In this study, every effort was made to minimise the above sources of error. Future research work should examine if it is possible to develop more accurate methods for the tracing of vertebral contour and measurement of intervertebral movements. This would be particularly helpful when it is desirable to determine intervertebral translations but the inverse kinematic method would not be reliable enough to predict such movements. The Hough transform (HT) is a powerful tool in computer vision for extracting objects in images by gathering evidence in a parameter space (LEAVERS, 1992). This method has been shown to be very promising for extracting the contours of vertebrae from radiographs (ZHENGet al., 2003).

5 Conclusions This was the first study that examined the feasibility of using a robotic engineering method to study the kinematics of the intervertebral movements of the lumbar spine. The inverse kinematic algorithm developed was found to be valid for predicting the rotational movements of the intervertebral joint for a given forward bending movement of the trunk. Such a prediction would require knowledge of the total flexion movement of the lumbar spine and the positions of the spinous processes, which could be determined by surface measurements. it is suggested that the technique could be used clinically to predict intervertebral rotations when it is clinically undesirable to measure such movements with radiographs or other methods carrying health risks. The technique will also be of value in biomechanical modelling when kinematic data are required. However, the study demonstrated that the method was unable reliably to predict intervertebral translation. The technique is of little value in the assessment of some spinal pathologies, such as spondylolysis or spondylolisthesis, when clinicians are required to determine if there is excessive translation of the intervertebral joint. Prediction of the translational movement is technically 745

difficult, as the magnitude o f the movement is small. Finally, it is concluded that furore research should further explore the use o f robotic theories for examining the kinematic mechanisms o f the lumbar spine. Acknowledgments" The authors wish to thank the Duchess o f Kent Children's Hospital, Hong Kong, for providing the radiographs. They also acknowledge the financial support o f the internal Competitive Research Grant (ICRG A-PD88) o f the Hong Kong Polytechnic University.

References ALLARD, R, STOKES, I. A. E, and BLANCHI, J. R (1995): 'Three-dimensional analysis of Human Movement' in 'Human Kinetics' (Champaign, IL, 1995) BIRYUKOVA,E. V, ROBY-BRAMI,A., FROLOV,A. A., and MOKHTARI, M. (2000): 'Kinematics of human arm reconstructed from spatial tracking system recordings', J. Biomech., 33, pp. 985-995 BROWN, R. H., BURSTEIN, A. H., NASH, C. L., and SCHOCK,C. C. (1976): 'Spinal analysis using a three-dimensional radiographic technique', J. Biomech., 9, pp. 355-365 CHAFFIN, D. B., ANDERSSON,G. B. J., and MARTIN, B. J. (1999): 'Occupational biomechanics', 3rd edn (Wiley InterScience, New York, 1999) CHEN, K., GIBLIN, R, and IRVING,A. (1999): 'Mathematical explorations with Matlab' (Cambridge University Press, 1999) CRAIG, J. J. (1989): 'Introduction to robotics: mechanics and control' (Addison-Wesley, 1989) DVORAK,J., PANJABI,M. M., CHANG,D. G., THEILER,K., and GROB, D. (1991): 'Functional radiographic diagnosis of the lumbar spine. Flexion-extension and lateral bending', Spine, 16, pp. 562-571 FUJIE, H., MABUCHI, K., WOO, S. L., LIVESAY,G. A., ARAI, S., and TSUKAMOTO, Y. (1993): 'The use of robotics technology to study human joint kinematics: A new methodology' J. Biomech. Eng., 115, pp. 211-217 GOLUB, G. H., and VAN LOAN, C. E (1996): 'The singular value decomposition and unitary matrixes', in 'Matrix computations, 3rd edn', (Johns Hopkins University Press, Baltimore, MD, 1996.) JUNG,E. S., CHOE,J., and KIN, S. H. (1994): 'PsychophysicaJ cost function of joint movement for arm reach posture prediction'. Proc 38th AnnuaJ Meeting Humaxl Factors and Ergonomics Society, 1, pp. 636-640 KAIGLE, A. M., POPE, M. H., FLEMING, B. C., and HANSSON, T. (1992): 'Technical note: a method for the intravital measurement of interspinous kinematics' J Biomeeh., 25, pp. 451-456 KANAYAMA,M., ABUMI, K., KANEDA,K., TADANO,S., and UKAI, T. (1996): 'Phase lag of the intersegmental motion in flexion-extension of the lumbar and lumbosacral spine: am in vivo study', Spine, 21, pp. 1416-1422 LEAVERS, V E (1992): 'Shape detection in computer vision using Hough transform' (Springer Verlag, 1992) LEE, R. Y. W. (1995): 'The biomechanical basis of spinal manual therapy'. PhD thesis, University of Strathclyde, Glasgow, UK LEE, R. Y. W., and EVANS, J. H. (1997): 'An in vivo study of the intervertebral movements produced by posteroanterior mobilisation', Clin. Biomech., 12, pp. 400-408 LEE, R. Y. W., and MUNN, J. (2000): 'Passive moment about the hip in straight leg raising', Clin. Biomech., 12, pp. 330-334 LEE, R. Y. W. (2001): 'Kinematics of rotational mobilisation of the lumbar spine', Clin. Biomech., 16, pp. 481-488 McCARTHY, J. M. (1990): 'Introduction to theoretical kinematics', (MIT Press, Cambridge, MA, 1990) McGILL, S. M., and NORMAN, R. W. (1986): 'Partitioning of the L4L5 dynamic moment into disc, ligamentous, and muscular components during lifting', Spine, 11, pp. 666-678 PANJABI,M. M., CHANG, D., and DVORAK, J. (1992a): 'An analysis of errors in kinematic paxaxneters associated with in vivo functional radiographs', Spine, 17, pp. 200-205 PANJABI, M. M., GOEL, V., OXLAND, T., TAKATA,K., DURANCEAU,J., KRAG, M., and PRICE, M. (1992b): 'Human lumbar vertebrae. Quantitative three-dimensional anatomy', Spine, 17, pp. 29%306 PEARCY, M. J. (1985): 'Stereo radiography of lumbar spine motion', Acta. Orthop. Scand., 56, pp. 1-45.

746

PEARCY, M. J., and HINDLE, R. J. (1989): 'New method for the non-invasive three-dimensional measurement of human back movement', Clin. Biomech., 4, pp. 73-79. PORTER, R. W., WICKS, M., and OTTEWELL,D. (1978): 'Measurement of the spinal canal by diagnostic ultrasound', J. Bone Joint Surg., 60B, pp. 481-484 ROSEN, M., BREEN, A., HAMANN,W., HARBER, P., JAYSON,M. I. V, KELLY, E., et al. (1994): 'Management guidelines for back pain'. Report of a clinical standards advisory group committee on back pain, (HMSO London, 1994), pp. 65-71 SCHAFFER,W. O., SPRATT,K. E, WEINSTEIN,J., LEHMANN,T. R., and GOEL, V (1990): 'The consistency and accuracy of roentgenograms for measuring sagittal translation in the lumbar vertebral motion segment. An experimental model', Spine, 15, 741-750 SOMMERIll, H. J., and MILLER,N. R. (1980): 'A technique for kinematic modeling of anatomical joints', J Biomeeh. Eng., 102, pp. 311-317 STOKES, I. A. E, BEVINS, T. M., and LUNN, R. A. (1987): 'Back surface curvature and measurement of lumbar spinal motion', Spine, 12, pp. 355-361 SUZUKI, S., YAMAMURO,T., SHIKATA,J., SHIMIZU, K., and IIDA, H. (1989): 'Ultrasound measurement of vertebral rotation in idiopathic scoliosis', d Bone Joint Surg., 71B, pp. 252-255 TAKAYANAGI, K., TAKAHASHI, K., YAMAGATA, M., MORIYA, H., KITAHARA, H., and TAMAKI, T. (2001): 'Using cineradiography for continuous dynamic-motion analysis of the lumbar spine', Spine, 26, pp. 1858-1865 TENCER, A. E, and MAYER, T. G. (1983): 'Soft tissue strain and facet face interaction in the lumbar intervertebral joint. Part I: Input data and computational technique', J. Biomech. Eng., 105, pp. 201-209 TILG, B., FISCHER, G., MODRE, R., HANSER, E, MESSNARZ, B., SCHOCKE, M., KREMSER, C., BERGER, T., HINTRINGER, E, and ROITHINGER, E X. (2002a): 'Model-based imaging of cardiac electrical excitation in humans' kinematics of human arm reconstructed from spatial tracking system recordings', IEEE Trans. Ned. Imag., 21, pp. 1031-1039. VAN DEN BOGERT, A. J., SMITH, G. D., and NIGG, B. M. (1994): 'In vivo determination of the anatomical axes of the ankle joint complex: am optimization approach', J. Biomech., 27, pp. 1477-1488 WANG, X. (1999): 'A behavior-based inverse kinematics algorithm to predict arm prehension postures for computer-aided ergonomic evaluation', J. Biomech., 32, pp. 453-460 WHITE, A. A., and PANJABI,M. M. (1978): 'Clinical biomechanics of the spine', (J.B. Lippincott, Philadelphia, 1978) YAMAMOTO,I., PANJABI,M. M., CRISCO, J. J., and OXLAND,T. (1989): 'Three-dimensional movements of the whole lumbar spine lumbosacral joint', Spine, 14, pp. 1256-1260. ZATSIORSKY,V. M. (1998): 'Kinematics of Human Motion' in 'Human Kinetics' (Champaign, 1998) ZHANG, X., and CHAFFIN, D. B. (1996): 'An optimization-based differential inverse kinematics approach for modeling three-dimensional dynamic postures during seated reaching movements', in 'Technical report TR96-15' (University of Michigan, MI, USA) ZHENG, Y., ALLEN, R., and NbXON, M. S. (2003): 'Lumbar spine visualisation based on kinematic analysis from videofluoroscopic imaging', Ned. Eng. Phys., 25, pp. 171-179 ZHOU, S. H., McCARTHY, I. D., McGREGOR, A. H., COOMBS, R. R., and HUGHES, S. R (2000): 'Geometrical dimensions of the lower lumbar vertebrae--analysis of data from digitised CT images', Eu~ Spine. J., 9, pp. 242-248

Author Biography RAYMOND LEE obtained his PhD from the Bioengineering Unit, The University Strathclyde, UK. He has held academic appointments in HonE Kong and Australia, and is currently am Associate Professor in the Department of Rehabilitation Sciences, The HonE Kong Polytechnic University, HonE Kong. He had previously been awarded the British Council Fellowship in Rehabilitation Engineering and the Development Fellowship by the Association of Commonwealth Universities. His research focuses on kinematics of the lumbar spine, movement coordination, the effects of back pain on spinal kinematics and the biomechanical basis of therapeutic procedures.

Medical & Biological Engineering & Computing 2004, Vol. 42