Real-time Deformation of Soft Tissue Internal Structure with Surface Profile Variations using Particle System

TL;DR

This paper presents a real-time method for simulating internal soft tissue deformation based on surface profile variations using a particle system, enabling efficient and mesh-free internal structure modeling during surgery.

Contribution

The novel approach models internal tissue deformation in real-time with a particle system that does not require 3D mesh generation, reducing computational complexity.

Findings

Handles up to 10,000 particles in 0.3 seconds

Reduces computational burden with simplified calculations

Does not require generating 3D meshes

Abstract

Intraoperative observation of tissue internal structure is often difficult. Hence, real-time soft tissue deformation is essential for the localization of tumor and other internal structures. We propose a method to simulate the internal structural deformations in a soft tissue with surface profile variations. The deformation simulation utilizes virtual physical particles that receive interaction forces from the surface and other particles and adjust their positions accordingly. The proposed method involves two stages. In the initialization stage, the three-dimensional internal structure of the surface mesh is uniformly sampled using the particle expansion and attracting-repelling force models whilst simultaneously building the internal particle connections. In the simulation stage, under surface profile variations, we simulate the internal structural deformation based on a deformation force model that uses the internal particle connections. The main advantage of this method is that it greatly reduces the computational burden as it only involves simplified calculations and also does not require generating three-dimensional meshes. Preliminary experimental results show that the proposed method can handle up to 10,000 particles in 0.3s.

Equations

$$\mathrm{d}{p}_i={\mathbf{F}}_i/m$$

$${\mathbf{F}}_i=\mu {\mathbf{F}}_{surface\to i}+\sum _{j\in {\Omega }_i}{\mathbf{F}}_{j\to i}$$

$${{\bf{F}}_{surface\rightarrow i}}=\left\{{\begin{array}[]{*{20}{c}}0&{if{\kern 1.0pt}{\kern 1.0pt}{{\bf{n}}_{i}}\cdot\;{{\bf{d}}_{surface\rightarrow i}}\leq 0}\\ {-{{\bf{n}}_{i}}}&{if{\kern 1.0pt}{\kern 1.0pt}{{\bf{n}}_{i}}\cdot\;{{\bf{d}}_{surface\rightarrow i}}>0}\end{array}}\right.$$

$$h=1/\left(\frac{\alpha d_{i,j}}{r_i+r_j}+(1-\alpha )\right)$$

$${\mathbf{F}}_{j\to i}=\left(h^6-h^3\right){\mathbf{d}}_{j\to i}$$

$$\sum _{j\in {\Omega }_i}\left|{\mathbf{F}}_{j\to i}\right|=\left|\sum _{j\in {\Omega }_i}{\mathbf{F}}_{j\to i}\right|+\left|{\mathbf{F}}_{compress,i}\right|$$

$$\begin{array}[]{l}\left|{{\bf{F}}_{compress,i}^{ideal}}\right|=\frac{{12}}{{{{\left({\frac{{0.5\times 0.75({r_{i}}+{r_{j}})}}{{{r_{i}}+{r_{j}}}}+0.5}\right)}^{6}}}}-\frac{{12}}{{{{\left({\frac{{0.5\times 0.75({r_{i}}+{r_{j}})}}{{{r_{i}}+{r_{j}}}}+0.5}\right)}^{3}}}}\;={\rm{8}}{\rm{.8257}}\end{array}$$

$$\mathrm{d}{r}_i=-\ln (1.0+\eta (\left|{\mathbf{F}}_{compress,i}\right|-\left|{\mathbf{F}}_{compress,i}^{ideal}\right|)$$

$${\Omega }_i=\{j|d_{i,j}<0.6(r_i+r_j)\}.$$

$${{\bf{A}}_{relax,i}}={\left[{\begin{array}[]{*{20}{c}}{{{({p_{j1}}-{p_{i}})}^{T}}}\\ \vdots\\ {{{({p_{jn}}-{p_{i}})}^{T}}}\end{array}}\right]_{n\times 3}},{{\bf{A}}_{current,i}}={\left[{\begin{array}[]{*{20}{c}}{{{({q_{j1}}-{q_{i}})}^{T}}}\\ \vdots\\ {{{({q_{jn}}-{q_{i}})}^{T}}}\end{array}}\right]_{n\times 3}}$$

$$\Vert {\mathbf{A}}_{relax,i}{\mathbf{R}}_i-{\mathbf{A}}_{current,i}\Vert$$

$$[\mathbf{U},\mathbf{S},{\mathbf{V}}^T]=\mathrm{svd}({\mathbf{A}}_{relax,i}^T{\mathbf{A}}_{current,i})$$

$${\mathbf{R}}_i=\mathbf{U}{\mathbf{V}}^T$$

$${\mathbf{F}}_{i\to j}=\sigma \left({\mathbf{R}}_i^T(p_j-p_i)-(q_j-q_i)\right)$$

Full text

1 Introduction

In the operating room, real-time intraoperative localization of the tumor and other internal structures of a tissue is essential and challenging. Compared to the surgical preparation and planning stage, when the structure of a diseased tissue can be accurately obtained by the combination of many diagnostics technologies, the intraoperative localization of internal structures is more difficult. For example, only 50% of non-palpable tumors are visible by intraoperative ultrasound in the breast [1]. In contrast, in the process of many surgeries, the tissue surface can be directly observed and be modeled by construction methods. For example, Conley et al used a laser scanner to reconstruct the skin surface of the breast [2]. Hence, it is important to obtain the internal structural deformations in a soft tissue according to surface profile variations.

Deformation simulation is a key problem in the computer graphics field and many approaches have been proposed [3][4]. Among these methods, particle system-based methods have been proven to be efficient and are widely used [5][6]. The main idea of particle system is to utilize virtual physical particles which interact with each other and move accordingly. Hence, particle system has an advantage of being more close to the reality.

For performing particle system-based calculation, it is a prerequisite to know which particles have direct connections between them. Many deformation methods use polygon meshes to represent the connections [7]. However, for a three-dimensional (3D) structure, the Delaunay-based tetrahedral meshes generation is quite time consuming [8][9]. Besides, the Delaunay-based methods need vertices as a prerequisite for performing mesh generation. Another type of deformation methods defines the connections by Euclidean distance thresholding [10], which greatly reduced the computation burden. For both types of deformation methods, the initial uniformly distributed particles with connections would improve the accuracy.

With the built connections, the mass spring damper system [11] is often employed for calculating the interaction forces between particles. However, the mass spring damper system has a significant drawback of not being able to take into account the variation in rotation angles directly. Although this drawback can be compensated by the implicit triangular stability in most cases, the mass spring damper system cannot handle situations with a large deformation.

In this paper, we employ this concept of particle system and utilize particles to represent physical units of a tissue. The main contribution of this paper are twofold: first, a novel sampling algorithm is proposed to construct virtual physical internal connections of a tissue with uniformly distributed particles, second; a novel local deformation description algorithm is proposed for generating deformation forces between particles. Both algorithms have been shown to be fast and efficient.

2 Uniformly Sampling

We want to down-sample the object internal structure with a given number of particles for reducing the computational burden of the subsequent deformation estimation. For this purpose, we choose uniformly sampling which can make the deformation simulation smooth. Further, because in a real world situation, physical particles only have direct interactions when they come in contact, and deformations are propagated through connections defined by this contacting relationships. Hence, it is necessary to build the virtual physical connections between particles for the deformation simulation.

Our uniformly sampling algorithm utilizes the conception of particle system, where particles receive forces from the object surface and their neighboring particles and move accordingly. For particle $p_{i}$, the relationship between its motion and the received forces is given as follows:

[TABLE]

[TABLE]

where ${\bf{F}}_{i}$ is sum of all forces it receives, ${\bf{F}}_{surface\rightarrow i}$ is the interaction force applied from the surface to particle $p_{i}$, ${\bf{F}}_{j\rightarrow i}$ is the interaction force applied from $p_{j}$ to $p_{j}$ ($i\neq j$), $\Omega_{i}$ is the set of neighboring particles of $p_{i}$, $m$ is the particle mass which is the same for all particles, and $\mu$ is a coefficient which controls the relative weights.

In the following of this section we will discuss the details of surface force calculation and how particles interact with each other.

2.1 Surface Force

The surface force ${\bf{F}}_{surface\rightarrow i}$ is introduced to keep particles inside of the object surface. If a particle is outside then the surface force will attract it to the object.

[TABLE]

where ${\bf{d}}_{surface\rightarrow i}$ is the unit vector from to the nearest surface point to ${{p}}_{i}$. ${\bf{n}}_{i}$ is the unit normal direction of the nearest surface point, we assume ${\bf{n}}_{i}$ points from inside to outside. As in Eq. (2), the weight coefficient $\mu$ is a large number to make sure the interaction forces between particles cannot push particles outside of the surface.

2.2 Particles Interaction

For uniformly sampling, we employed the concept that some entities repel or attract each other to avoid creating a density that is a too high or too low. We created a simple physical model by considering each particle as a balloon that is inflated gradually. The balloons will repel each other when they come in contact. In this model we need to solve two problems: (1) how to calculate the interaction force. and (2) how to inflate or deflate the balloons.

The interaction force applied from particle $p_{j}$ to $p_{j}$ is based on the simplified intermolecular forces (IMFs) model.

[TABLE]

[TABLE]

where ${{\bf{d}}_{j\rightarrow i}}$ is the unit vector from $p_{j}$ to $p_{i}$, $d_{i,j}$ is the distance between $p_{i}$ and $p_{j}$, $r_{i}$ and $r_{j}$ are the imaginary balloon radii, $\alpha$ is a coefficient. We choose $\alpha=0.5$.

The simplified IMFs model (5) makes interaction force zero when $d_{i,j}=r_{i}+r_{j}$. When $d_{i,j}<r_{i}+r_{j}$ , ${\bf{F}}_{j\rightarrow i}\cdot{{\bf{d}}_{j\to i}}>0$ indicates particles will repel each other. When $d_{i,j}>r_{i}+r_{j}$ , ${\bf{F}}_{j\rightarrow i}\cdot{{\bf{d}}_{j\to i}}<0$ indicates particles will attract each other. When $p_{i}$ is far away from $p_{j}$, their interaction force will gradually decrease to zero.

In our algorithm, each balloon $p_{i}$ gradually dilates by imaginary inflation, but the particle interaction forces may compress the balloon. The particle interaction force applied on $p_{i}$ can be divided into two parts (see (6)). The first part on the right side of (6) is involved in the motion of $p_{i}$, and the second part compresses the balloon.

[TABLE]

It is not appropriate to make the end point of the convergence progress near $d_{i,j}=r_{i}+r_{j}$, where the interaction forces are near zero so particles take much more time to reach the goal positions. Hence, in our design, the algorithm converges at $d_{i,j}\approx 0.75(r_{i}+r_{j})$, where particles still repelling each other.

Considering a close-packing of equal spheres, which is a ideal uniform sampling pattern and each sphere has 12 adjacent spheres. Presuming that each balloon stops changing its radius at this ideal situation, hence:

[TABLE]

Consequently, $r_{i}$ updates according to the compression force:

[TABLE]

where $\eta$ is a coefficient suggesting the hardness of each balloon.

According to (8), when $\left|{{{\bf{F}}_{compress,i}}}\right|<\left|{{\bf{F}}_{compress,i}^{ideal}}\right|$, the radius $r_{i}$ becomes larger. According to (5), $p_{i}$ will have larger interaction forces with neighboring points than before, which tends to make $\left|{{{\bf{F}}_{compress,i}}}\right|$ larger. Similarly, when $\left|{{{\bf{F}}_{compress,i}}}\right|>\left|{{\bf{F}}_{compress,i}^{ideal}}\right|$, our algorithm tends to make $\left|{{{\bf{F}}_{compress,i}}}\right|$ smaller. This forms a negative feedback system, with the result that all variables in our algorithm, including the points distribution, become stable.

For particle $p_{i}$, we assume it is connected to each particles in its neighboring points set

[TABLE]

Fig. 1 shows the uniformly sampling process. Our method generates particles on both surface and inside. Those surface particles will be used as the control points in the subsequent deformation simulation and their position changes according to the change of surface.

3 Deformation Simulation

Our deformation simulation algorithm also utilizes particle system but with a different interaction force model. We consider the initial particles distribution and connections obtained by our uniformly sampling method as the relax state, which means there is no interaction force between particles. If surface profile varies, the surface particles, which are used as the control points, will change their positions accordingly and apply forces to internal particles. The surface particles deformation is propagated deeper into the internal structure and the final new distribution should minimize the energy of the particle system. To date, the mass spring damper system is wildly used for deformation simulation.

However, mass spring damper system has a main problem of not being able to take into account angles variation directly, as shown in Fig. 2(a). To overcome this drawback, we presented a method to calculate interaction forces by taking into account all deformation information. The idea of our method is straight forward, as shown in Fig. 2(b). For a particle $p_{i}$ and its connected particles $p_{j},j\in{\Omega_{i}}$, the change of their relative positions may cause interaction force. It is worth noting that the relative position variations are affected by both rigid rotation and shape deformation, we need to first get rid of the effect caused by rigid rotation.

Denoting particles at the relax state as $p$ and those at the current state as $q$, the related local shape at particle $i$ can be described in matrix form:

[TABLE]

We find a rotation matrix ${\bf{R}}_{i}$, which minimizes

[TABLE]

Hence, ${\bf{R}}_{i}$ can be found as follows:

[TABLE]

[TABLE]

The deformation force applied from particle $i$ to $j$ is:

[TABLE]

where $\sigma$ is a coefficient suggesting the stiffness of the local tissue.

After calculations of all points we can get the sum of forces that a particle receives, and then we move particles according to the motion model (1).

4 Experiments

The algorithms proposed in this paper were implemented in C++ and VTK running on a 2.60 GHz Intel Core i7 processer. Our experiments first analyzed the algorithms runtime, and then evaluated the accuracy based on real breast tumor data.

4.1 Runtime

Both algorithms are based on particle system with linear time complexity. We evaluated their runtime with different number of particles varying from 1,000 to 10,000.

As shown in Fig. 3(a), the uniformly sampling algorithm took less than 3 seconds for processing up to 10,000 particles. Because the uniformly sampling algorithm is only needed at the initialization stage, the real-time performance relies mostly on the deformation algorithm. As shown in Fig. 3(b), one iteration of the deformation algorithm took an average of less than 60ms for handling up to 10,000 particles, which is very fast. It should be noted that the required number of iterations depends on how large the relative deformation is between two states, hence it is not appropriate to give an exact total runtime needed for deformation simulation. In our experiments we found that for consistently tracking a deforming object in realtime, it usually needed less than five iterations to change from one state to a new stable state. This suggests that for 10,000 particles, the average deformation runtime is around 0.3s.

4.2 Experiments on Breast Tumor Data

Breast tumor MR images are used to evaluate the accuracy of our algorithms. We calculated the effect of breast skin deformation on the motion of the breast tumor. We segmented and modelled the chest wall, breast skin and tumor by using 3D Slicer. The breasts skin deformation were obtained by aligning intraoperative prone and supine breast MR imaging, as shown in Fig.4(a). The deformation process are shown in Fig.4(b)-(d). The tumor was considered as a rigid object, whose center was represented with a particle and its orientation was adjusted according to the connected neighboring particles. The accuracy comparison results are shown in Fig. 5. Two patients’ data were used and the error distance of the center of the deformed prone tumor and supine tumor were 2.8 mm and 3.5 mm respectively.

5 Conclusion

In this paper, we present a particle system-based method for real-time deformation of soft tissues, which achieves high accuracy with an error of several millimeters. Because the calculation of each particle is independent, the algorithms is high parallel. Future work will involve the implementation a GPU parallel computing version of the algorithms. Future work should also involve improving the particle initialization algorithm to allow generating higher density particles near region of interest to achieve high accuracy.