A new ChainMail approach for real-time soft tissue simulation

ABSTRACT

This paper presents a new ChainMail method for real-time soft tissue simulation. This method enables the use of different material properties for chain elements to accommodate various materials. Based on the ChainMail bounding region, a new time-saving scheme is developed to improve computational efficiency for isotropic materials. The proposed method also conserves volume and strain energy. Experimental results demonstrate that the proposed ChainMail method can not only accommodate isotropic, anisotropic and heterogeneous materials but also model incompressibility and relaxation behaviors of soft tissues. Further, the proposed method can achieve real-time computational performance.

KEYWORDS:

• ChainMail method
• real-time performance
• surgery simulation
• soft tissue deformation

Introduction

Soft tissue simulation is a challenging research topic in the field of surgery simulation. Surgery simulation requires realistic and real-time modeling of tool-tissue interactions; however, it is difficult to meet both of these conflict requirements.¹ The literature on soft tissue deformation can be divided into 2 classes, with one focused on real-time capability such as mass-spring model (MSM) and the other focused on accurate deformation such as finite element method (FEM).^{2, 3} The former is computationally efficient but does not allow accurate modeling of soft tissue material properties, whereas the latter is computationally extensive and formulations must be simplified to reduce runtime computation.⁴ Various techniques have been reported to improve the computational performance of FEM. The matrix condensation reduces the computational time by confining the full computation of a volume mesh to the surface nodes; however, this simplification significantly reduces the simulation accuracy.⁵ The pre-computation reduces the computational time based on a set of pre-computed spatial derivatives; however, it does not allow any changes on model topology.⁶ The explicit total Lagrangian FEM reduces the computational load by eliminating iteratively solving a large system of equations; however, the solution is only stable under strictly small time steps due to the use of explicit time integration.^6-8 The GPU (Graphics Processing Unit) acceleration facilitates the computational performance of FEM; however, it relies on hardware and does not solve computational problem fundamentally.^9,10 In general, with most of the existing techniques in FEM, the improved computational performance is achieved by sacrificing the accuracy of FEM modeling.

ChainMail is a modeling method alternative to the above methods. The basic concept of this method is a chain element (mass point) enforces a bounding region for each of its neighboring chain elements to control their movements. The first ChainMail method, the 3D ChainMail, was proposed by Gibson for real-time soft tissue simulation.¹¹ This method has the significant advantage in computation and can handle topology change. It is also simple in implementation and stable in numerical iteration.¹² Although improvements on the ChainMail method have been studied, the method still lacks the capability of modeling soft tissues' incompressibility and relaxation behaviors due to the lack of volume and strain energy conservation.^13,14

This paper presents a new ChainMail method for real-time soft tissue simulation. The proposed ChainMail represents an alternative to the FEM for real-time soft tissue simulation. It can handle isotropic, anisotropic and heterogeneous materials by simply changing the material parameters of each chain element. A new time-saving scheme is developed for isotropic materials to improve the computational efficiency. The proposed ChainMail method also achieves the conservation of volume and strain energy for modeling the incompressibility and relaxation behaviors of soft tissues. Experiments and comparison analysis have been conducted to evaluate the performance of the proposed method.

Traditional ChainMail method

From the viewpoint of continuum mechanics, the ChainMail method is equivalent to a spring system. As shown in Figure 1, a spring of length ${\displaystyle l_0}$ at the rest state can be compressed to a minimum length of ${\displaystyle l_{min}}$ and extended to a maximum length of ${\displaystyle l_{max}}$. Therefore, the movement of the spring is bounded within the minimum compression limit ${\displaystyle l_{min}}$ and the maximum extension limit ${\displaystyle l_{max}}$.

Figure 1. The spring of rest length ${\displaystyle l_0}$ with its minimum compression length ${\displaystyle l_{min}}$ and maximum extension length ${\displaystyle l_{max}}$.

Traditional ChainMail method, such as the 3D ChainMail and the Generalized ChainMail, employs ChainMail bounding regions (CBRs) to regulate the movements of chain elements in the object.^11,13As illustrated in Figure 2, chain element ${\displaystyle P_j\left(x_{\text{j}},y_j,z_{\text{j}}\right)}$ is moved to ${\displaystyle P_j^{*}(x_{\text{j,i}}^{*},y_{\text{j,i}}^{*},z_{\text{j,i}}^{*})}$ while ${\displaystyle P_i\left(x_{\text{i}},y_{\text{i}},z_{\text{i}}\right)}$ is moved to ${\displaystyle P_i^{*}\left(x_{\text{i}}^{*},y_{\text{i}}^{*},z_{\text{i}}^{*}\right)}$. The ${\displaystyle {\text{CBR}}_{\text{j,i}}}$ for chain element ${\displaystyle P_j}$ with respect to (w.r.t) the new position ${\displaystyle P_i^{*}}$ of chain element ${\displaystyle P_i}$ in the Generalized ChainMail is given by Eq. (1(1) ${\text{CBR}}_{\text{j,i}}=\begin{array}{l}\{x_{\text{j,i min}}\le x_{\text{j,i}}^{*}\le x_{\text{j,i max}};y_{\text{j,imin}}\le y_{\text{j,i}}^{*}\\ \le y_{\text{j,i max}};z_{\text{j,i min}}\le z_{\text{j,i}}^{*}\le z_{\text{j,i max}\}}\end{array}$(1) ).

Figure 2. The ${\displaystyle {\text{CBR}}_{\text{j,i}}}$, shown as the dashed bounding box, for chain element ${\displaystyle P_j}$ w.r.t chain element ${\displaystyle P_i^{*}}$.

(1) ${\text{CBR}}_{\text{j,i}}=\begin{array}{l}\{x_{\text{j,i min}}\le x_{\text{j,i}}^{*}\le x_{\text{j,i max}};y_{\text{j,imin}}\le y_{\text{j,i}}^{*}\\ \le y_{\text{j,i max}};z_{\text{j,i min}}\le z_{\text{j,i}}^{*}\le z_{\text{j,i max}\}}\end{array}$(1) where ${\displaystyle x_{\text{j,i min}}}$ and ${\displaystyle x_{\text{j,i max}}}$ are the minimum compression and maximum extension limits in the x axis, i.e.(2) ${\displaystyle x_{\text{j,i min}}=x_{\text{i}}^{*}+\left({\alpha }_{\text{min}}\Delta x_{\text{j,i}}-\beta \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right);x_{\text{j,i max}}=x_{\text{i}}^{*}+\left({\alpha }_{\text{max}}\Delta x_{\text{j,i}}+\beta \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)}$(2) where ${\displaystyle {\alpha }_{\text{min}}}$,${\displaystyle {\alpha }_{\text{max}}}$ and ${\displaystyle \beta }$ are material parameters, and ${\displaystyle \Delta x_{\text{j,i}}}$, ${\displaystyle \Delta y_{\text{j,i}}}$ and ${\displaystyle \Delta z_{\text{j,i}}}$ are geometric distances between chain elements ${\displaystyle P_i}$ and ${\displaystyle P_j}$ w.r.t the x, y and z axes at the rest state, i.e., ${\displaystyle \Delta x_{\text{j,i}}=\left|x_{\text{j}}-x_{\text{i}}\right|}$. The limits in the other 2 directions (y and z) are expressed in a similar manner. As the 3 material parameters are constant throughout the entire volumetric object, the traditional ChainMail cannot model non-linear material properties such as anisotropy and inhomogeneity.

Proposed TSVE-ChainMail

A new ChainMail method named the TSVE-ChainMail, time-saving volume-energy conserved ChainMail, is presented. The TSVE-ChainMail can handle various material properties, improve computational efficiency for isotropic materials, and conserve objects' volume and strain energy.

New ChainMail bounding region

The limits ${\displaystyle x_{\text{j,i min}}}$ and ${\displaystyle x_{\text{j,i max}}}$ of the new ${\displaystyle {\text{CBR}}_{\text{j,i}}}$ in the proposed TSVE-ChainMail are defined as(3) ${\displaystyle \begin{array}{ll}{x}_{\text{j,i min}}=x_{\text{i}}^{*}+\left(\Delta x_{\text{j,i}}-\frac{{\alpha }_{\text{i}}+{\alpha }_{\text{j}}}{2}\left(\Delta x_{\text{j,i}}+\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)& x_{\text{j,i max}}=x_{\text{i}}^{*}+\left(\Delta x_{\text{j,i}}+\frac{{\alpha }_{\text{i}}+{\alpha }_{\text{j}}}{2}\left(\Delta x_{\text{j,i}}+\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)\end{array}}$(3) where ${\displaystyle {\alpha }_{\text{i}}}$ and ${\displaystyle {\alpha }_{\text{j}}}$ are material parameters of chain elements ${\displaystyle P_i}$ and ${\displaystyle P_j}$, and they are corresponded to the spring stiffness ${\displaystyle k}$. For isotropic materials (${\displaystyle {\alpha }_{\text{i}}={\alpha }_{\text{j}}=\alpha )}$, Eq. (3(3) ${\displaystyle \begin{array}{ll}{x}_{\text{j,i min}}=x_{\text{i}}^{*}+\left(\Delta x_{\text{j,i}}-\frac{{\alpha }_{\text{i}}+{\alpha }_{\text{j}}}{2}\left(\Delta x_{\text{j,i}}+\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)& x_{\text{j,i max}}=x_{\text{i}}^{*}+\left(\Delta x_{\text{j,i}}+\frac{{\alpha }_{\text{i}}+{\alpha }_{\text{j}}}{2}\left(\Delta x_{\text{j,i}}+\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)\end{array}}$(3) ) may be simplified to(4) ${\displaystyle \begin{array}{ll}{x}_{\text{j,i min}}=x_{\text{i}}^{*}+\left(\left(1-\alpha \right)\Delta x_{\text{j,i}}-\alpha \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)& x_{\text{j,i max}}=x_{\text{i}}^{*}+\left(\left(1+\alpha \right)\Delta x_{\text{j,i}}+\alpha \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)\right)\end{array}}$(4)

Anisotropic and heterogeneous materials can be modeled by setting different parameters in different directions and regions respectively, which was unable to achieve with the traditional ChainMail.

Time saving scheme for isotropic materials

A time saving scheme (TSS) is developed for isotropic materials of constant material parameter ${\displaystyle \alpha }$, where each chain element is considered only once in every iteration, leading to improved computational efficiency. Consider the 3 general cases of chain elements' movement illustrated in Fig. 3: Case (1) Following ${\displaystyle P_i}$'s movement to ${\displaystyle P_i^{*}}$, ${\displaystyle P_j}$ and ${\displaystyle P_k}$ are moved to ${\displaystyle P_j^{*}}$ and ${\displaystyle P_k^{*}}$ respectively (Fig. 3(a)). Then${\displaystyle P_j^{*}}$ and ${\displaystyle P_k^{*}}$ stay within their respective CBRs; Case (2) Following ${\displaystyle P_k}$'s movement to ${\displaystyle P_k^{*}}$, ${\displaystyle P_l}$ is moved to ${\displaystyle P_l^{*}}$ (Fig. 3(b)). Then ${\displaystyle P_j^{*}}$ and ${\displaystyle P_l^{*}}$ stay within their respective CBRs; Case (3) In addition to ${\displaystyle P_i}$, ${\displaystyle P_j}$, ${\displaystyle P_k}$ and ${\displaystyle P_l}$'s movements, ${\displaystyle P_n}$ is moved to ${\displaystyle P_n^{*}}$ following ${\displaystyle P_i}$'s movement to ${\displaystyle P_i^{*}}$, whereas ${\displaystyle P_m}$ is moved to ${\displaystyle P_m^{*}}$ following ${\displaystyle P_n}$'s movement (Fig. 3(c)). Then ${\displaystyle P_l^{*}}$ and ${\displaystyle P_m^{*}}$ stay within their respective CBRs.

Figure 3. The chain element ${\displaystyle P_i}$ is moved to ${\displaystyle P_i^{*}}$ (blue arrow) while others follow its movement; the rest state is shown in black and the new positions are shown in red; the solid line indicates that the 2 elements are connected and one is moved directly w.r.t the other, whereas the dotted line indicates that the elements are connected but one is not moved w.r.t the other.

[Proof of Case (1)] Case (1) will be verified if ${\displaystyle P_j^{*}(x_{\text{j,i}}^{*},y_{\text{j,i}}^{*},z_{\text{j,i}}^{*})}$ stays within the ${\displaystyle {\text{CBR}}_{\text{j,k}}}$ w.r.t ${\displaystyle P_k^{*}}$, and ${\displaystyle P_k^{*}}$ stays within the ${\displaystyle {\text{CBR}}_{\text{k,j}}}$ w.r.t ${\displaystyle P_j^{*}}$. Since the exact positions of ${\displaystyle P_j^{*}}$ and ${\displaystyle P_k^{*}}$ are unknown, the maximum and minimum limits (2 limits), ${\displaystyle P_{k,max}^{*}\left(x_{\text{k,i max}},y_{\text{k,i max}},z_{\text{k,i max}}\right)}$ and ${\displaystyle P_{k,min}^{*}\left(x_{\text{k,i min}},y_{\text{k,i min}},z_{\text{k,i min}}\right)}$, are used for the verification of Case (1). Denote the ${\displaystyle {\text{CBR}}_{\text{j,k}}}$ at the 2 limits of ${\displaystyle P_k^{*}}$ by ${\displaystyle \overline{{\text{CBR}}_{\text{j,k}}}}$. Thus,(5) ${\displaystyle \overline{{\text{CBR}}_{\text{j,k}}}=\left\{\overline{x_{\text{j,k min}}}\le x_{\text{j,i}}^{*}\le \overline{x_{\text{j,k max}}};\overline{y_{\text{j,k min}}}\le y_{\text{j,i}}^{*}\le \overline{y_{\text{j,k max}}};\overline{z_{\text{j,k min}}}\le z_{\text{j,i}}^{*}\le \overline{z_{\text{j,k max}}}\right\}}$(5)

Substituting ${\displaystyle x_{\text{i}}^{*}=x_{\text{k,i min}}-\left(1-\alpha \right)\Delta x_{\text{k,i}}+\alpha \left(\Delta y_{\text{k,i}}+\Delta z_{\text{k,i}}\right)}$ into the lower limit ${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{j,i min}}=x_{\text{i}}^{*}+\left(1-\alpha \right)\Delta x_{\text{j,i}}-\alpha \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)}$ yields(6) ${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{k,i min}}-\left(1-\alpha \right)\Delta x_{\text{k,i}}+\alpha \left(\Delta y_{\text{k,i}}+\Delta z_{\text{k,i}}\right)+\left(1-\alpha \right)\Delta x_{\text{j,i}}-\alpha \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)}$(6)

Using the geometric distance ${\displaystyle \Delta x_{\text{j,i}}=\left|x_{\text{j}}-x_{\text{i}}\right|}$, Eq. (6(6) ${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{k,i min}}-\left(1-\alpha \right)\Delta x_{\text{k,i}}+\alpha \left(\Delta y_{\text{k,i}}+\Delta z_{\text{k,i}}\right)+\left(1-\alpha \right)\Delta x_{\text{j,i}}-\alpha \left(\Delta y_{\text{j,i}}+\Delta z_{\text{j,i}}\right)}$(6) ) may be simplified to(7) ${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{k,i min}}+\left(1-\alpha \right)\Delta x_{\text{j,k}}-\alpha \left(\Delta y_{\text{j,k}}+\Delta z_{\text{j,k}}\right)=\overline{x_{\text{j,k min}}}}$(7)

A similar calculation can be made for ${\displaystyle \overline{x_{\text{j,k max}}}}$, ${\displaystyle \overline{y_{\text{j,k min}}}}$,${\displaystyle \overline{y_{\text{j,k max}}},\overline{z_{\text{j,k min}}},\overline{z_{\text{j,k max}}}}$, and for ${\displaystyle {\text{CBR}}_{\text{k,j}}}$. Therefore, it is demonstrated that ${\displaystyle P_j^{*}}$ stays within the ${\displaystyle {\text{CBR}}_{\text{j,k}}}$ w.r.t ${\displaystyle P_k^{*}}$, and vice versa.

[Proof of Case (2)] Case (2) will be verified if ${\displaystyle P_l^{*}(x_{\text{l,k}}^{*},y_{\text{l,k}}^{*},z_{\text{l,k}}^{*})}$ stays within the ${\displaystyle {\text{CBR}}_{\text{l,j}}}$ w.r.t ${\displaystyle P_j^{*}}$, and ${\displaystyle P_j^{*}}$ stays within the ${\displaystyle {\text{CBR}}_{\text{j,l}}}$ w.r.t ${\displaystyle P_l^{*}}$. Denote the ${\displaystyle {\text{CBR}}_{\text{l,j}}}$ at the 2 limits of ${\displaystyle P_j^{*}}$ by ${\displaystyle \overline{{\text{CBR}}_{\text{l,j}}}}$. Thus,(8) $\overline{{\text{CBR}}_{\text{l,j}}}=\begin{array}{l}\{\overline{x_{\text{l,j min}}}\le x_{\text{l,k}}^{*}\le \overline{x_{\text{l,j max}}};\overline{y_{\text{l,j min}}}\le y_{\text{l,k}}^{*}\le \overline{y_{\text{l,j max}}};\overline{z_{\text{l,j min}}}\le z_{\text{l,k}}^{*}\le \overline{z_{\text{l,j max}}}\}\end{array}$(8)

Considering the lower limit ${\displaystyle x_{\text{l,k}}^{*}\ge x_{\text{l,k min}}=x_{\text{k,i}}^{*}+\left(1-\alpha \right)\Delta x_{\text{l,k}}-\alpha (\Delta y_{\text{l,k}}+\Delta z_{\text{l,k}}).}$ Substituting ${\displaystyle x_{\text{k,i}}^{*}}$ by ${\displaystyle x_{\text{k,i min}}=\overline{x_{\text{j,k min}}}-\left(1-\alpha \right) \Delta x_{\text{j,k}}+\alpha \left(\Delta y_{\text{j,k}}+\Delta z_{\text{j,k}}\right)}$ from Eq. (7(7) ${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{k,i min}}+\left(1-\alpha \right)\Delta x_{\text{j,k}}-\alpha \left(\Delta y_{\text{j,k}}+\Delta z_{\text{j,k}}\right)=\overline{x_{\text{j,k min}}}}$(7) ) yields(9) ${\displaystyle x_{\text{l,k}}^{*}\ge \overline{x_{\text{j,k min}}}-\left(1-\alpha \right)\Delta x_{\text{j,k}}+\alpha \left(\Delta y_{\text{j,k}}+\Delta z_{\text{j,k}}\right)+\left(1-\alpha \right)\Delta x_{\text{l,k}}-\alpha \left(\Delta y_{\text{l,k}}+\Delta z_{\text{l,k}}\right)}$(9)

Using the geometric distance, Eq. (9(9) ${\displaystyle x_{\text{l,k}}^{*}\ge \overline{x_{\text{j,k min}}}-\left(1-\alpha \right)\Delta x_{\text{j,k}}+\alpha \left(\Delta y_{\text{j,k}}+\Delta z_{\text{j,k}}\right)+\left(1-\alpha \right)\Delta x_{\text{l,k}}-\alpha \left(\Delta y_{\text{l,k}}+\Delta z_{\text{l,k}}\right)}$(9) ) may be simplified to(10) ${\displaystyle x_{\text{l,k}}^{*}\ge \overline{x_{\text{j,k min}}}+\left(1-\alpha \right)\Delta x_{\text{l,j}}-\alpha \left(\Delta y_{\text{l,j}}+\Delta z_{\text{l,j}}\right)}$(10)

Based on Eqs. (1(1) ${\text{CBR}}_{\text{j,i}}=\begin{array}{l}\{x_{\text{j,i min}}\le x_{\text{j,i}}^{*}\le x_{\text{j,i max}};y_{\text{j,imin}}\le y_{\text{j,i}}^{*}\\ \le y_{\text{j,i max}};z_{\text{j,i min}}\le z_{\text{j,i}}^{*}\le z_{\text{j,i max}\}}\end{array}$(1) ) and (7(7) ${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{k,i min}}+\left(1-\alpha \right)\Delta x_{\text{j,k}}-\alpha \left(\Delta y_{\text{j,k}}+\Delta z_{\text{j,k}}\right)=\overline{x_{\text{j,k min}}}}$(7) ),${\displaystyle x_{\text{j,i}}^{*}\ge x_{\text{j,i min}}}$ and ${\displaystyle x_{\text{j,i}}^{*}\ge \overline{x_{\text{j,k min}}}}$. At the lower limit of ${\displaystyle x_{\text{j,i}}^{*}}$, ${\displaystyle \overline{x_{\text{j,k min}}}=x_{\text{j,i min}}}$,(11) ${\displaystyle x_{\text{l,k}}^{*}\ge x_{\text{j,i min}}+\left(1-\alpha \right)\Delta x_{\text{l,j}}-\alpha \left(\Delta y_{\text{l,j}}+\Delta z_{\text{l,j}}\right)=\overline{x_{\text{l,j min}}}}$(11)

A similar calculation can be made for the other boundary limits, and for ${\displaystyle {\text{CBR}}_{\text{j,l}}}$. Therefore, it is demonstrated that ${\displaystyle P_l^{*}}$ stays within the ${\displaystyle {\text{CBR}}_{\text{l,j}}}$ w.r.t ${\displaystyle P_j^{*}}$, and vice versa.

[Proof of Case (3)] Case (3) will be verified if ${\displaystyle P_m^{*}(x_{\text{m,n}}^{*},y_{\text{m,n}}^{*},z_{\text{m,n}}^{*})}$ stays within the ${\displaystyle {\text{CBR}}_{\text{m,l}}}$ w.r.t ${\displaystyle P_l^{*}}$, and ${\displaystyle P_l^{*}}$ stays within the ${\displaystyle {\text{CBR}}_{\text{l,m}}}$ w.r.t ${\displaystyle P_m^{*}}$. Denote the ${\displaystyle {\text{CBR}}_{\text{l,m}}}$ at the 2 limits of ${\displaystyle P_m^{*}}$ by ${\displaystyle \overline{{\text{CBR}}_{\text{l,m}}}}$. Thus,(12) $\overline{{\text{CBR}}_{\text{l,m}}}=\begin{array}{l}\{\overline{x_{\text{l,m min}}}\le x_{\text{l,k}}^{*}\le \overline{x_{\text{l,m max}}};\overline{y_{\text{l,m min}}}\le y_{\text{l,k}}^{*}\le \overline{y_{\text{l,m max}}};\overline{z_{\text{l,m min}}}\le z_{\text{l,k}}^{*}\le \overline{z_{\text{l,m max}}}\}\end{array}$(12)

Substituting ${\displaystyle x_{\text{j,i min}}=\overline{x_{\text{m,j} \text{min}}}-\left(1-\alpha \right)\Delta x_{\text{m,j}}+\alpha \left(\Delta y_{\text{m,j}}+\Delta z_{\text{m,j}}\right)}$ in the lower limit given by Eq. (11(11) ${\displaystyle x_{\text{l,k}}^{*}\ge x_{\text{j,i min}}+\left(1-\alpha \right)\Delta x_{\text{l,j}}-\alpha \left(\Delta y_{\text{l,j}}+\Delta z_{\text{l,j}}\right)=\overline{x_{\text{l,j min}}}}$(11) )(13) ${\displaystyle x_{\text{l,k}}^{*}\ge \overline{x_{\text{m,j} \text{min}}}-\left(1-\alpha \right)\Delta x_{\text{m,j}}+\alpha \left(\Delta y_{\text{m,j}}+\Delta z_{\text{m,j}}\right)+\left(1-\alpha \right)\Delta x_{\text{l,j}}-\alpha \left(\Delta y_{\text{l,j}}+\Delta z_{\text{l,j}}\right)}$(13)

Similar to ${\displaystyle x_{\text{l,k}}^{*}\ge \overline{x_{\text{l,j min}}}}$ in Eq. (11(11) ${\displaystyle x_{\text{l,k}}^{*}\ge x_{\text{j,i min}}+\left(1-\alpha \right)\Delta x_{\text{l,j}}-\alpha \left(\Delta y_{\text{l,j}}+\Delta z_{\text{l,j}}\right)=\overline{x_{\text{l,j min}}}}$(11) ),${\displaystyle x_{\text{m,n}}^{*}\ge \overline{x_{\text{m,j min}}}}$. At lower limit of ${\displaystyle x_{\text{m,n}}^{*}}$, ${\displaystyle \overline{x_{\text{m,j min}}}=x_{\text{m,n min}}}$. Thus,(14) ${\displaystyle x_{\text{l,k}}^{*}\ge x_{\text{m,n min}}+\left(1-\alpha \right)\Delta x_{\text{l,m}}-\alpha \left(\Delta y_{\text{l,m}}+\Delta z_{\text{l,m}}\right)=\overline{x_{\text{l,m min}}}}$(14)

A similar calculation can be made for the other boundary limits, and for ${\displaystyle {\text{CBR}}_{\text{m,l}}}$. Therefore, it is demonstrated that ${\displaystyle P_l^{*}}$ stays within the ${\displaystyle {\text{CBR}}_{\text{l,m}}}$ w.r.t ${\displaystyle P_m^{*}}$, and vice versa. From the above, we can draw the following 3 remarks. For any 3 neighboring elements A, B and C,

[Remark 1] if A and B are moved w.r.t C, then A and B stay within their respective CBRs (${\displaystyle P_i^{*}}$, ${\displaystyle P_j^{*}}$ and ${\displaystyle P_k^{*}}$ in Fig. 3(a)).

[Remark 2] if only A is moved w.r.t C while B and C have been previously moved as per Remark 1, then A and B stay within their respective CBRs (${\displaystyle P_j^{*}}$, ${\displaystyle P_k^{*}}$ and ${\displaystyle P_l^{*}}$ in Fig. 3(b)).

[Remark 3] if A and C as well as B and C have been previously moved as per Remark 2, then A and B stay within their respective CBRs (${\displaystyle P_j^{*}}$, ${\displaystyle P_l^{*}}$ and ${\displaystyle P_m^{*}}$ in Fig. 3(c)).

The three remarks demonstrate that the chain elements always satisfy the CBRs w.r.t each other in isotropic materials. Hence, each chain element can be considered only once at each iteration, saving the computational time. It should be noted that boundary conditions can be enforced by specifying displacement values to chain elements at the boundary of the calculation domain.

Volume and strain energy conservation

Since most biological soft tissues are incompressible and there is a change in strain energy when a soft object is deformed by an external force, the conservation of volume and strain energy are introduced for accurate soft tissue simulation by adding an additional position adjustment ${\displaystyle \Delta P_i}$.¹⁵ This position adjustment is derived based on the conditions of conservation of volume and energy.^16,17(15) ${\displaystyle \Delta P_i=\frac{w_{i}C\left(P_0,\dots ,P_n\right)}{{{\displaystyle \sum }}_{j=1}^4{w}_j{\left|{\nabla }_{P_j}C\left(P_0,\dots ,P_n\right)\right|}^2}{\nabla }_{P_i}C\left(P_0,\dots ,P_n\right)}$(15) (16) ${\displaystyle C\left(P_0,\dots ,P_n\right)=C_{volume}\left(P_0,P_1,P_2,P_3\right)=\frac{1}{6}\left|\left(P_1-P_0\right)·\left(\left(P_2-P_0\right)\times \left(P_3-P_0\right)\right)\right|-V_0}$(16) (17) ${\displaystyle C\left(P_0,\dots ,P_n\right)=C_{energy}\left(P_0,P_1,P_2,P_3\right)=\frac{1}{2}{\displaystyle \sum }_{i=0}^3{k}_i{\overrightarrow{({\mathbf{\text{P}}}_{\mathbf{\text{i}}}{\mathbf{\text{P}}}_{\mathbf{\text{C}}}}-l_i)}^2}$(17) where ${\displaystyle P_0,P_1,P_2,P_3}$ are the vertices of a tetrahedron,${\displaystyle w_i}$ and ${\displaystyle w_j}$ are the inverse of masses of respective chain elements, ${\displaystyle V_0}$ is the volume of the tetrahedron at rest state, ${\displaystyle P_C}$ is the barycenter of the tetrahedron, ${\displaystyle \overrightarrow{{\mathbf{\text{P}}}_{\mathbf{\text{i}}}{\mathbf{\text{P}}}_{\mathbf{\text{C}}}}}$ and ${\displaystyle l_i}$ are the current and the rest length of ${\displaystyle \overrightarrow{{\mathbf{\text{P}}}_{\mathbf{\text{i}}}{\mathbf{\text{P}}}_{\mathbf{\text{C}}}}}$, and ${\displaystyle k_i}$ is the stiffness of chain element ${\displaystyle P_i}$.

Implementation results and evaluations

A prototype surgery simulation system has been implemented with the proposed TSVE-ChainMail. Experiments have been conducted to evaluate the performance of the proposed method in terms of soft tissue material properties, computational time, volume conservation and strain energy conservation.

Fig. 4 illustrates the deformation of isotropic, anisotropic and heterogeneous materials. The cubic model was deformed evenly for an isotropic material in Fig. 4(a), deformed more significantly in one direction than the others in Fig. 4(b), and deformed differently at different regions in Fig. 4(c).

Figure 4. Isotropic, anisotropic and heterogeneous deformations.

The computational performance was evaluated on an Intel(R) Core(TM) i7-4700 [email protected] PC. Experiments were conducted under the same conditions to compare the timing performances with and without the use of TSS for isotropic materials. As illustrated in Fig. 5, the computational time with TSS was less than that without TSS for isotropic materials. A computational gain of 29.4% has been achieved for the real-time visual refresh rate (30Hz) in terms of the number of moving chain elements. The volume conservation has been achieved with the proposed ChainMail. As shown in Fig. 6, the proposed ChainMail conserved the object's volume considerably better than the traditional ChainMail, which had a significant volume loss after deformation. The volume change after deformation was 0.71% for the TSVE-ChainMail whereas it was 14.31% for the traditional ChainMail.

Figure 5. Computational time of the ChainMails with and without TSS for isotropic materials.

Figure 6. Comparison between the ChainMails in terms of volume conservation.

Interactive deformation of virtual human organs with force feedback has been achieved with the proposed method. Fig. 7(a) illustrates the deformation process of a volumetric kidney model with a virtual probe. The strain energy conservation is presented in Fig. 7(b)–(d): the kidney model returned to its original state with the proposed ChainMail whereas the traditional ChainMail failed to do so. It was also noticed that a significant visual improvement has been achieved with the proposed method.

Figure 7. Comparison between the ChainMails in terms of strain energy conservation: the kidney model with the TSVE-ChainMail can return to its original state once the external force is released whereas the traditional ChainMail fails to do so.

Conclusion

This paper presents a time-saving volume-energy conserved ChainMail method for real-time soft tissue simulation. This method allows different materials to be assigned to different chain elements to handle various material behaviors. A time-saving scheme is developed to improve computational efficiency for isotropic materials, and volume and strain energy conservation are proposed for realistic soft tissue deformation. Results demonstrate that the proposed method cannot only handle isotropic, anisotropic and heterogeneous materials but also model soft tissues' incompressibility and relaxation behaviors. The achieved performance is outperformed than that of the traditional ChainMail. Future research work will focus on extending the proposed method to model complex surgical operations.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.