Vascular deformation modeling scheme
Efficient vascular deformation modeling scheme plays an important role in virtual reality based vascular interventional surgery training system since it determines the visual deformation result and real-time performance directly. In this section, we first discuss a hybrid geometric blood vessel model which processes collision detection and deformation calculation separately, and then introduce a position-based vascular deformation method with volume conservation.
Hybrid geometric blood vessel model
A hybrid geometric blood vessel model which consists of a relative simplified triangular mesh and a complicated tetrahedral mesh is employed in this paper to improve the efficiency of collision detection between catheter and vascular wall. Figure 2a depicts the conceptual diagram of the proposed hybrid blood vessel model where the tetrahedral mesh represents vascular wall with specific thickness inspired by the vascular anatomy model, as shown in Fig. 2b. The red tetrahedral mesh is mainly used for accurate deformation calculation, while the green triangular mesh represents the ghost surface (never rendering) that tightly attached to the inner vascular wall. This triangular mesh is mainly designed for collision detection between the catheter and the blood vessel and further preventing the catheter from penetrating the blood vessel. Compared with previous related work [9], employing position-based dynamics technique with volume conservation and hash table based spatial acceleration scheme for real-time vessel deformation simulation is our main concern in this paper. Therefore, we simply model the catheter with tetrahedral mesh (yellow cylinder) mainly for further efficient collision detection and convenient experimental results validation.
The motion states of a flexible catheter and its contact information with vascular wall are also important in this paper. In our implementation, the catheter is also modeled with tetrahedral mesh and its tip has four degrees of freedom (DOF) just as a real catheter, including three degrees of freedom in rotation (i.e. rotation around x, y and z axis) and one translation along the catheter axis direction. The rest parts are controlled by elastic potential energy produced by the collision response (collision between the catheter and the ghost triangular mesh). In this way, a catheter can move freely in the blood vessel with continuous transmission commands governed by an experienced surgeon.
Position-based vascular deformation with volume conservation
Compared with the widely-used mass-spring system, position-based dynamic approach is unconditional stable even solved with a large time step. The original studies are mainly used to model the deformation of triangular mesh like cloth simulation [19, 20]. Cloth can be constructed with a large amount of triangles easily and its constraints mainly consist of stretching and bending constraints. As we mentioned above, the elements in our deformation calculation are tetrahedrons rather than triangles so that volume conservation constraint should be also taken into account as well as stretching and bending constraints to prevent the vascular wall from the issue of penetration and makes the deformation much more realistic. However, the volume conservation constraint is just mentioned briefly without detailed derivation in the original study [18]. To implement the position-based vascular deformation with volume conservation, a detailed derivation of the constraint should be given firstly. Let p be the vertices concatenation \({\bigl [p_1^T, p_2^T, p_3^T, p_4^T\bigr ]}^{T}\) of a tetrahedron. Given a vertex \(p_i\) with mass \(m_i\), where \(i \in (1,2,3,4)\) is the index number of a vertex, as well as j in the following equation. Then, the displacement of the vertex \(p_i\) by projection can be represented as
$$\begin{aligned} \Delta p_i =-\frac{C(p)}{\sum _{j=1}^4 w_j {|\nabla _{p_j} C(p)|}^{2}} w_i \nabla _{p_i} C(p) \end{aligned}$$
(1)
where \(w_i=1 / m_i\), C(p) and \(\nabla _p C(p)\) mean the constraint function and its gradient respectively. According to the reference [19], the volume conservation constraint of a single tetrahedron with rest volume \(V_0\) can be described as
$$\begin{aligned} C(p_1, p_2, p_3, p_4)=\frac{1}{6}\left( (p_2 - p_1) \times (p_3 - p_1)\right) \cdot (p_4 - p_1)- V_0 \end{aligned}$$
(2)
And our ultimate goal is to deduce every unknown factor in Eq. (1). The volume of a tetrahedron can be calculated with
$$\begin{aligned} V(p_1, p_2, p_3, p_4)&= \frac{1}{6}\left( (p_2 - p_1) \times (p_3 - p_1)\right) \cdot (p_4 - p_1) \nonumber \\&= \frac{1}{6}\left( (p_2 \times p_3 - p_2 \times p_1 - p_1 \times p_3 + p_1 \times p_1)\right) \cdot (p_4 - p_1) \nonumber \\&= \frac{1}{6}\left( (p_2 \times p_3)\cdot p_4 - (p_2 \times p_1)\cdot p_4 - (p_1 \times p_3)\cdot p_4- (p_2 \times p_3)\cdot p_1\right) \end{aligned}$$
(3)
According to Eq. (1), we have to calculate \(\nabla _{p_1}V,\nabla _{p_2}V,\nabla _{p_3}V\) and \(\nabla _{p_4}V\) respectively. In order to make the derivation clearly, the basic gradient principles of a combination of cross and dot product can be represented as follows, where \(a,b,c \in R^3\)
$$\begin{aligned} \nabla _c(a \times b)\cdot c=\, & {} \nabla _c\left( (a_2b_3-a_3b_2)c_1 + (a_3b_1-a_1b_3)c_2 +(a_1b_2-a_2b_1)c_3\right) \nonumber \\=\, & {} a \times b,\end{aligned}$$
(4)
$$\begin{aligned} \nabla _c(a \times c)\cdot b=\, & {} \nabla _c\left( (a_2c_3-a_3c_2)b_1 + (a_3c_1-a_1c_3)b_2 +(a_1c_2-a_2c_1)b_3\right) \nonumber \\=\, & {} b \times a,\end{aligned}$$
(5)
$$\begin{aligned} \nabla _c(c \times a)\cdot b=\, & {} \nabla _c\left( (c_2a_3-c_3a_2)b_1 + (c_3a_1-c_1a_3)b_2 +(c_1a_2-c_2a_1)b_3\right) \nonumber \\=\, & {} a \times b. \end{aligned}$$
(6)
With the principles (4)–(6), we can deduce the corresponding gradient of Eq. (3) easily
$$\begin{aligned} \nabla _{p_1}V=\, & {} \frac{1}{6}(p_2 \times p_4 - p_3 \times p_4 - p_2 \times p_3) \nonumber \\=\, & {} \frac{1}{6}(p_2 \times p_4 - p_3 \times p_4 - p_2 \times p_3 - p_2 \times p_2) \nonumber \\=\, & {} \frac{1}{6}((p_4 - p_2)\times (p_3 - p_2)) \end{aligned}$$
(7)
Similarly, the gradient of other corners can be obtained
$$\begin{aligned} \nabla _{p_2}V=\, & {} \frac{1}{6}\left( (p_3-p_1)\times (p_4-p_1)\right) \\ \nabla _{p_3}V=\, & {} \frac{1}{6}\left( (p_4-p_1)\times (p_2-p_1)\right) \\ \nabla _{p_4}V=\, & {} \frac{1}{6}\left( (p_2-p_1)\times (p_3-p_1)\right) \end{aligned}$$
Then according to the above equations and Eqs. (4)–(7), the denominator in Eq. (1) can be represented as
$$\begin{aligned} \sum ^4_{j=1}w_j|\nabla p_j C(p)|^2 &= w_1\left| \frac{1}{6}(p_4-p_2)\times (p_3 -p_2)\right| ^2 + w_2\left| \frac{1}{6}(p_3-p_1)\times (p_4 -p_1)\right| ^2 \nonumber \\&\quad + w_3\left| \frac{1}{6}(p_4-p_1)\times (p_2 -p_1)\right| ^2 + w_4\left| \frac{1}{6}(p_2-p_1)\times (p_3 -p_1)\right| ^2 \end{aligned}$$
(8)
Now we have all the necessary factors in Eq. (1), which also means the desired volume conservation constraint function is achieved. With the volume conservation constraint mentioned above, realistic deformation result of a surgery training system can be improved greatly and our detailed experimental results will be given in the section of results and discussion.
Hash table based spatial acceleration
The blood vessel model in a surgery training system usually has a large amount of tetrahedral elements, which make the vascular deformation calculation with volume conservation constraint much more time-consuming and may further limit the real-time performance of the system seriously. To guarantee a high real-time performance without losing deformation accuracy, several advanced algorithms and efficient data structures can be used. Inspired by the previous related work [21], a hash table based spatial acceleration algorithm is employed in our study to reduce computation overhead so as to significantly improve the real-time performance of our surgery training system. However, the original approach is used to process rapid collision detection between soft bodies, while in this paper we mainly use the hash value to quickly determine the vascular deformation region and then perform partial calculation which can improve real-time performance greatly.
Calculating vascular deformation partially is the core idea of the spatial acceleration approach. Figure 3 depicts the spatial acceleration approach with a small piece of blood vessel model, in which the wireframe and the filled box represent the inactive and active bounding boxes respectively. In the initialization stage, all tetrahedral elements of a blood vessel model are first classified into several adjacent child bounding boxes whose states are inactive by default and each accompanies with a unique hash value. The hash value is extremely important since it determines which bounding box encloses the catheter tip. As we mentioned above, a catheter can move freely in the blood vessel under the control of a skilled surgeon, thus a hash value of the catheter tip need to be calculated in every time-step. A bounding box will be activated when its hash value matches the hash value of the catheter tip, namely the active bounding box encloses the catheter tip. And further, all tetrahedral elements in the active bounding box are also active and will be used for deformation calculation, which means that an active bounding box is a basic calculation unit in our implementation.
During the simulation, some particles are fixed to prevent the blood vessel from falling under the gravity and moving freely under the pulled force. In our study, the particles which are too close to the top or bottom of the active bounding box are fixed. However, the catheter moves from one bounding box to another dynamically and when the catheter tip just enters a new bounding box with small displacement, the vascular deformation can hardly occur because of the fixed particles. To address this issue and guarantee more smooth transition when the catheter tip moves from one bounding box to another, we expanded the child bounding box a little bigger than the regular one when classifying the tetrahedral elements into a specific child bounding box. As a result, a tetrahedral element may be existed in several adjacent bounding boxes at the same time because of the shared space. The larger the shared space is, the smoother the transition results can be obtained while the more computation overhead is required. The detailed hash table based spatial acceleration approach will be discussed below.
The entire simulation domain is enclosed by an axis-aligned bounding box (AABB) with \(g^{max}\) and \(g^{min}\) denoting its corresponding minimum and maximum coordinates, then the entire bounding box should be partitioned into several regular child bounding boxes with size d. Each bounding box accompanies with a hash value which stored in an array, the length of the array is \(N_{array} = N_x \cdot N_y \cdot N_z\), and
$$\begin{aligned} (N_x, N_y, N_z) = \biggl (\biggl \lceil \frac{g^{max}_x-g^{min}_x}{d} \biggr \rceil , \biggl \lceil \frac{g^{max}_y-g^{min}_y}{d} \biggr \rceil , \biggl \lceil \frac{g^{max}_z-g^{min}_z}{d} \biggr \rceil \biggr ) \end{aligned}$$
(9)
where \(N_x, N_y\) and \(N_z\) denote the number of bounding boxes in the three perpendicular axes respectively. Once given a specific bounding box with three axis-based indices (i, j, k), we can use the following function T (i, j, k) to uniquely map the three-dimensional bounding box into an element of the continuous one-dimensional array.
$$\begin{aligned} T(i,j,k) = i + j\times N_x + k\times N_x \cdot N_y, (0,0,0) \le (i,j,k) \le (N_x,N_y,N_z) \end{aligned}$$
(10)
In function T (i, j, k), a triple loop is used to map the bounding box into the array one by one continuously without overlapping. And the regular child bounding box can be calculated with
$$\begin{aligned} \left\{ \begin{array}{lll} e^{min}(i,j,k) &{}=&{} g^{min}+d\times (i,j,k) \\ e^{max}(i,j,k) &{}=&{} e^{min}(i,j,k)+(d,d,d) \end{array} \right. \end{aligned}$$
(11)
According to the reference [21], the corresponding hash value of a bounding box can be determined with
$$\begin{aligned} h(T(i,j,k))=\bigl [(i\times u)\oplus (j\times v)\oplus (k\times w)\bigr ] \ mod \ N_{array} \end{aligned}$$
(12)
Here u, v and w are three large prime numbers to ensure a unique hash value for each bounding box.
Since a regular child bounding box is designed to be a basic computation unit, we have to determine all the tetrahedral elements that a regular child bounding box contains. Each regular child bounding box is expanded with the size \(span_x=span_y=span_z=d/3\) to guarantee a much more smooth transition during the simulation. The elongation span is an experienced value and it can balance computation efficiency and visual deformation result quite well.
During the simulation, the catheter moves in the blood vessel model freely and the position of the catheter tip is \((x_t,y_t,z_t)\) at a given time t and then we can calculate the corresponding hash value with
$$\begin{aligned} h_{Tip} = \biggl [\biggl (\biggl \lfloor {\frac{x_t-g^{min}_x}{d}}\biggl \rfloor \times u\biggr ) \! \oplus \! \biggl (\biggl \lfloor {\frac{y_t-g^{min}_y}{d}}\biggl \rfloor \times v\biggr ) \! \oplus \! \biggl (\biggl \lfloor {\frac{z_t-g^{min}_z}{d}}\biggl \rfloor \times w\biggr )\biggr ] \ mod \ N_{array} \end{aligned}$$
(13)
The parameters u, v and w have the same meaning and value as in Eq. (12). Then the child bounding box with the same hash value of the catheter tip will be activated and used for precise collision detection and deformation calculation.