**6.2 Haptic rendering method for the rigid stage of the spine models**

In real-time haptic simulation, users can only interact with the spine model by manipulating a rigid virtual object considered as a probe on the computer screen. At present, a simple probe such as a sphere is used in this study. In the rigid stage, a common haptic rendering method is used for the aforementioned types of spine models. Since this interaction carries out at a high update rate of 1 kHz, the chosen haptic rendering method needs to be

Development of a Detailed Human Spine Model with Haptic Interface 183

Fig. 19. Complete haptic simulation process in the system

Fig. 20. Collision between the probe and a vertebra of the spine model

reasonable and effectively computational. As presented above, the haptic rendering method includes two phases: collision detection and collision response. For the collision detection, the simple algorithm proposed by James Arvo (1990) is utilized to check intersection between the probe and the spine model through axis-aligned bounding boxes (AABBs). Figure 20 illustrates all AABBs of a vertebra intersecting with the spherical probe.

For collision response, the algorithm developed by Gao and Gibson (2006) is used to determine force feedback. An important prerequisite in this step is that intersecting points with the spherical probe have to be specified. Figure 21 shows intersecting points during colliding process between the probe and a vertebra.

Fig. 18. Procedure of haptic rendering

reasonable and effectively computational. As presented above, the haptic rendering method includes two phases: collision detection and collision response. For the collision detection, the simple algorithm proposed by James Arvo (1990) is utilized to check intersection between the probe and the spine model through axis-aligned bounding boxes (AABBs).

For collision response, the algorithm developed by Gao and Gibson (2006) is used to determine force feedback. An important prerequisite in this step is that intersecting points with the spherical probe have to be specified. Figure 21 shows intersecting points during

Figure 20 illustrates all AABBs of a vertebra intersecting with the spherical probe.

colliding process between the probe and a vertebra.

Fig. 18. Procedure of haptic rendering

Fig. 19. Complete haptic simulation process in the system

Fig. 20. Collision between the probe and a vertebra of the spine model

Development of a Detailed Human Spine Model with Haptic Interface 185

In this stage, the haptic interface has two components: haptic rendering and an FE solver. The function of the haptic rendering component is to provide haptic feedback to the user when the user virtually deforms the spine model. The FE solver evaluates the deformation

To solve the FE deformation problem, an explicit Newmark integration method is employed. The theory of elasticity consists of a set of differential equations that describe the state of stress, strain and displacement of each point within an elastic deformable body. A relationship between the deformation of the object and the exerted forces can be established by synthesizing those equations. The equation governing the resulting deformation from

3 3 *n n* where u is the 3n-dimensional nodal displacement vector; and are the velocity and acceleration vectors respectively; F is the external force vector; M is the mass matrix; D is the damping matrix; R(u) represents the internal force vectors due to deformation; n is the

If deformation is assumed to be small, a linear elastic model can be used to approximate

*u x*

*u v y x*

 

where x, y and z are the independent variables of the Cartesian frame, and u, v and w are

Because the internal force vector is linear with respect to the nodal displacement vector, this leads to a constant stiffness matrix, which can be pre-computed. However, this solution does not consider geometric non-linearity and is only fit therefore for linear and small deformation situations. The linear strain model introduces distortions for large global deformations. In the model the deformation of the spine is rather large and therefore the linear strain approximation is not suitable. We model the deformation using the Green

22 2 1

*u uvw x xxx*

*u v uu vv ww y x xy xy xy*

 

*x*

*xy*

the corresponding displacement variables at the given point.

2 *<sup>x</sup>*

*xy*

*Mu Du R u F* ( ) (3)

(4)

(5)

(6)

(7)

**6.3.1 Finite element spine model** 

number of nodes in the FE model.

external forces is:

linear strain:

strain as follows:

of the FE spine model according to the user's force input.

Fig. 21. Intersecting points betweenthe probe and the vertebra

After these intersecting points are found, the final step in this phase is to calculate the reaction forces generated from the vertebrae of the spine model. Computing the reaction forces of the vertebrae is concerned with finding force magnitudes which the probe applied to the intersecting points on the contact surface areas of the vertebrae and force directions at those points. The magnitudes and directions of reaction forces can be computed by the equations as follows:

$$\begin{aligned} \left| F \right| = \frac{\sum\_{i=0}^{n} \left( k^\* \left. s\_i - d \right| \* \left. \dot{s}\_i \right)}{S\_{\text{MaxSec}}} \right. \tag{1} \end{aligned} \tag{1}$$

$$S\_{\text{MaxSec}}$$

where k is spring constant, s is the displacement of intersected point, d is the damping factor, is the velocity of the intersected point, is the maximum sectional area of the probe.

$$\vec{F} = \frac{\sum\_{i=0}^{n} \vec{N}\_i \, \ast \, s\_i}{n} \tag{2}$$

where n is the number of surface sampling points inside probe, Ni as the normal vector of surface sampling point and si as the depth.

#### **6.3 Haptic rendering method for the compliant stage of the spine models**

After the users locate a specific vertebra where he/she wishes to apply force, they can then press the PHANToM stylus button and push or drag the vertebra in any direction to make the whole spine model deform. Once the stylus button is pressed, the system switches to another haptic rendering algorithm that uses the stretched-spring model. A virtual spring is set up connecting the vertebra and the haptic cursor. The spring has two hook points: one is on the vertebra and the other is the haptic cursor. Both of the hook points displace during spine deformation. The force magnitude is determined by the length of the virtual spring and stiffness while the force direction depends on the vector of the virtual spring.

#### **6.3.1 Finite element spine model**

184 Haptics Rendering and Applications

After these intersecting points are found, the final step in this phase is to calculate the reaction forces generated from the vertebrae of the spine model. Computing the reaction forces of the vertebrae is concerned with finding force magnitudes which the probe applied to the intersecting points on the contact surface areas of the vertebrae and force directions at those points. The magnitudes and directions of reaction forces can be computed by the

*i i*

(1)

(2)

\* \*

*ks ds*

*MaxSec*

*S* 

*SMaxSec*

where k is spring constant, s is the displacement of intersected point, d is the damping factor, is the velocity of the intersected point, is the maximum sectional area of the probe.

0

where n is the number of surface sampling points inside probe, Ni as the normal vector of

After the users locate a specific vertebra where he/she wishes to apply force, they can then press the PHANToM stylus button and push or drag the vertebra in any direction to make the whole spine model deform. Once the stylus button is pressed, the system switches to another haptic rendering algorithm that uses the stretched-spring model. A virtual spring is set up connecting the vertebra and the haptic cursor. The spring has two hook points: one is on the vertebra and the other is the haptic cursor. Both of the hook points displace during spine deformation. The force magnitude is determined by the length of the virtual spring

*i*

 

*F*

**6.3 Haptic rendering method for the compliant stage of the spine models** 

and stiffness while the force direction depends on the vector of the virtual spring.

*n*

\*

*N s*

*n*

*i i*

0

*i*

*F*

*n*

Fig. 21. Intersecting points betweenthe probe and the vertebra

equations as follows:

surface sampling point and si as the depth.

In this stage, the haptic interface has two components: haptic rendering and an FE solver. The function of the haptic rendering component is to provide haptic feedback to the user when the user virtually deforms the spine model. The FE solver evaluates the deformation of the FE spine model according to the user's force input.

To solve the FE deformation problem, an explicit Newmark integration method is employed. The theory of elasticity consists of a set of differential equations that describe the state of stress, strain and displacement of each point within an elastic deformable body. A relationship between the deformation of the object and the exerted forces can be established by synthesizing those equations. The equation governing the resulting deformation from external forces is:

$$M\ddot{\boldsymbol{\mu}} + D\dot{\boldsymbol{\mu}} + R(\boldsymbol{\mu}) = \boldsymbol{F} \tag{3}$$

#### 3 3 *n n*

where u is the 3n-dimensional nodal displacement vector; and are the velocity and acceleration vectors respectively; F is the external force vector; M is the mass matrix; D is the damping matrix; R(u) represents the internal force vectors due to deformation; n is the number of nodes in the FE model.

If deformation is assumed to be small, a linear elastic model can be used to approximate linear strain:

$$
\varepsilon\_{\rm x} = \frac{\partial u}{\partial \mathbf{x}}\tag{4}
$$

$$\gamma\_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial \mathbf{x}} \tag{5}$$

where x, y and z are the independent variables of the Cartesian frame, and u, v and w are the corresponding displacement variables at the given point.

Because the internal force vector is linear with respect to the nodal displacement vector, this leads to a constant stiffness matrix, which can be pre-computed. However, this solution does not consider geometric non-linearity and is only fit therefore for linear and small deformation situations. The linear strain model introduces distortions for large global deformations. In the model the deformation of the spine is rather large and therefore the linear strain approximation is not suitable. We model the deformation using the Green strain as follows:

$$
\varepsilon\_{\mathbf{x}} = \frac{\partial u}{\partial \mathbf{x}} + \frac{\mathbf{1}}{\mathbf{2}} \left[ \left( \frac{\partial u}{\partial \mathbf{x}} \right)^2 + \left( \frac{\partial v}{\partial \mathbf{x}} \right)^2 + \left( \frac{\partial w}{\partial \mathbf{x}} \right)^2 \right] \tag{6}
$$

$$\mathcal{N}\_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} + \left[\frac{\partial u \partial u}{\partial x \partial y} + \frac{\partial v \partial v}{\partial x \partial y} + \frac{\partial w \partial w}{\partial x \partial y}\right] \tag{7}$$

Development of a Detailed Human Spine Model with Haptic Interface 187

relationships can be interpolated using the least-squares method and expressed in terms of polynomial functions. Figure 22 illustrates the graph of a dynamic property of vertebra T1 under forward forces in sagittal plane of the spine model. In these figures, each series of markers presents each dynamic property of one vertebra obtained in the simulations and the

**Translation of vertebra T1 under Forward External Force**

0 100 200 300 400 500 600

83 42 <sup>2</sup> *y eF eF eF* 6.3857 1.425 1.7442 0.09279 *zz z*

**Fz (N)**

Initially, without pressing the stylus button, the user can touch and explore the spine model since it is considered to be rigid throughout. Then, the user can control the probe to apply external forces in any arbitrary direction onto any vertebra by pressing the stylus button to observe the locomotion and study dynamic properties of the spine model. Figure 23 shows the haptic real-time simulation of the finite element spine model under lateral bending. During this real-time simulation, the displacement, rotation and deformation of the spine can be recorded and then input to the offline FEA application for further detailed

stress/strain and deformation behaviour analysis of IVDs as well as vertebrae.

**7.2 Studying stress/strain of the finite element spine model in a prone position** 

As presented above, all dynamic properties of the spine model recorded during the haptic simulation can be provided as input of offline FEA simulator for further detailed analysis. In


**7. Applications** 

Fig. 22. Relative translation y of T1 versus forward force

**7.1 Haptic simulation of the finite element spine model** 






**y(mm)**





0

continuous line closely fit to that series is the corresponding interpolated line.

The other 4 terms of the strain are defined similarly. The algorithm equations are:

$$
\dot{\mu}\_{n+1} = \dot{\mu}\_n + \frac{\Delta t}{2} (\ddot{\mu}\_n + \dddot{\mu}\_{n+1}) \tag{8}
$$

$$
\mu\_{n+1} = \mu\_n + \dot{\mu}\_n \Delta t + \frac{\ddot{\mu}\_n}{2} \Delta t^2 \tag{9}
$$

$$\dot{M}\ddot{u}\_{n+1} = F\_{n+1} - R(\mu\_{n+1}) - D\dot{u}\_{n+1} \tag{10}$$

Zhuang's mass concentration method is employed to approximate the distributed mass in real-time (Zhuang & Canny, 2000). To avoid inverting a large sparse matrix at each time step, Zhuang's method approximates the force matrix M with a diagonal matrix. Diagonalization results in approximating the mass continuum as concentrated masses at each nodal point of the finite element model. This simplifies the nonlinear system of equations used to form a set of independent algebraic equations that can be solved with relatively high updating rate.

A prototype system has been implemented. Although the FE model of the spine is a simplified beam model, the updating rate of FEM evaluation is still far below the 1000 Hz required for effective vibration-free haptic rendering. The FE model has 472 elements and 473 nodes and the time step of the explicit integration is 0.02 second. The updating rate is approximately 200 Hz using a PC graphic workstation. To avoid vibration due to the updating rate discrepancy, a force interpolation method (Zhuang & Canny, 2000) was adopted to smoothen the force feedback. Haptic feedback between two simulated states is linearly interpolated at a suitably high frequency. However, at time tn, we do not have the information of tn+1. Zhuang's solution uses an intentional delay of an FEM deformation time step, up to 1/100s. User experimentation found that such a small lag in time is within the tolerance of human perception for virtual interaction with soft objects.

#### **6.3.2 Multi-body spine model**

Different from the finite element spine model whose locomotion is directly calculated in the haptic simulation, the multi-body spine model in the compliant stage utilize the preinterpolated displacement-force functions of all vertebrae to conveniently compute the movement of the spine. Based on the magnitude of the haptic external forces applied by the users, dynamic properties of all vertebrae can be easily calculated via these functions and the locomotion of the whole spine model can be rapidly observed.

The purpose of interpolating displacement-force functions of all vertebrae is to facilitate computational process of dynamic properties of the spine model during haptic simulation. To do that, some constraints are imposed on the spine model. The pelvis is fixed in 3D space. Then, constant forces are applied on each specific vertebra in the thoracic region in each axisaligned direction during simulation. The force magnitude is gradually increased with an equal increment in subsequent simulations. Corresponding to each value of force, dynamic characteristics of all vertebrae (e.g. translation, rotation) can be automatically obtained using the plots in LifeMOD as a reference. The dynamic properties are recorded after the spine model is stabilized. Based on these recorded dynamic properties, displacement-force relationships can be interpolated using the least-squares method and expressed in terms of polynomial functions. Figure 22 illustrates the graph of a dynamic property of vertebra T1 under forward forces in sagittal plane of the spine model. In these figures, each series of markers presents each dynamic property of one vertebra obtained in the simulations and the continuous line closely fit to that series is the corresponding interpolated line.

Fig. 22. Relative translation y of T1 versus forward force
