Kalman Filter Finite Element Method for Real-Time Soft Tissue Modeling

Soft tissue modelling plays a significant role in surgery simulation as well as surgical procedure planning and training. However, it is a challenging research task to satisfy both physical realism and real-time simulation for soft tissue deformation. The finite element method (FEM) is a representative strategy for modelling of soft tissue deformation with highly physical realism. However, it suffers from expensive computations, unable to meet the requirement of real-time simulation. This paper proposes a novel method by combining FEM with the Kalman filter for real-time and accurate modelling of soft tissue deformation. The novelty of this method is that soft tissue deformation is formulated as a filtering identification process to online estimate soft tissue deformation from local measurement of displacement. To construct the discrete system state equation for filtering estimation, soft tissue deformation is discretised based on elastic theory in the space domain by FEM and is further discretised in the time domain by using the Wilson- $\theta $ implicit integration to solve the dynamic equilibrium equation of FEM deformation modelling. Subsequently, a Kalman filter is developed for online estimation and analysis of soft tissue deformation according to local measurement of displacement. Interactive tool-tissue interaction with haptic feedback is also achieved for surgery simulation. The presented method significantly improves the computational performance of the traditional FEM, but still maintains a similar level of accuracy. It not only achieves the real-time performance, but also exhibits the similar deformation behaviours as the traditional FEM and enables the use of large time steps to improve the simulation efficiency.


I. INTRODUCTION
Soft tissue modelling is a key issue in surgery simulation as well as surgical procedure planning and training. It requires soft tissues to respond to the external force in physical realism and real time. However, it is difficult to satisfy these two conflicting requirements [1], [2]. Currently, the existing methods for modelling of soft tissue deformation can be divided into two general classes [3], [4]. One focuses on real-time computational performance. The typical examples include the massspring [5] and chainmail [6], [7]. These methods are easy to implement and less time-consuming in calculation. However, they do not allow accurate modelling of material properties and the increase of the number of springs or chains also leads to a stiffer system. Recently, the traditional mass-spring The associate editor coordinating the review of this manuscript and approving it for publication was Pietro Savazzi . model was improved using flexion springs for real-time shape restoration [8]. However, this improved mass-spring model is based on surface representation without the interior geometry, unsuitable for complex surgical operations such as tearing and cutting. The other class focuses on accurate deformation modelling based on continuum mechanics. The typical example is the finite element method (FEM) [9], in which the mechanical behaviours of the soft tissue are accurately characterized based on rigorous strain-stress relationships by discretizing the soft tissue model into finite volumetric elements [10], [11]. However, FEM is computationally expensive, even for linear elasticity [12], [13].
Considerable research efforts have focused on improving the FEM computing performance. The technique of matrix condensation decreases the FEM computing time by limiting the entire calculation of volumetric elements to boundary surface nodes, but at the cost of sacrificing the accuracy [14]. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ The pre-computation method pre-calculates the equilibrium solutions to reduce the FEM computing time [15], which leads to a complex and time-consuming pre-computation process. The explicit FEM disentangles the sophisticated computing process of FEM by lumping masses and internal forces to individual nodes [14]. However, its solution is stable subject to small steps only. The GPU (Graphics Processing Unit) acceleration is a technique to facilitate the FEM computing performance via GPU [16], [17]. Nevertheless, this technique requires the hardware component of GPU for implementation. The total Lagrangian formulation combines the precomputation of spatial derivatives with the explicit integration to reduce the FEM computing time [18], [19]. However, it does not accept topological alternations in the model and its solution is stable subject to small steps only. The meshless method simplifies the FEM model discretization by discarding the mesh topology information [20], [21]. However, this method still requires the computation of node-to-node adjacency at each time step, causing an extra computational load. Further, it has difficulty in handling sparsely sampled regions and lacks theoretical error bounds on numerical integration [12]. The model order reduction techniques such as POD (Proper Orthogonal Decomposition) [22] and PGD (Proper Generalized Decomposition) [23] reduce the FEM computational load by approximating the complete displacement domain using a set of variables of much lower dimensions. However, this improvement is attained at the cost of reduced accuracy. FEM-based machine learning is also a method to achieve the real-time performance [24]. Although the increase of the sample size can increase the accuracy, it also increases the computational load for the learning process. Generally, most of the existing methods enhance the FEM computing performance at the cost of losing the modelling accuracy [20], [25]. The state-of-art survey on the existing modelling techniques of soft tissue deformation including various FEM improvements can be found in [12], [25], [26]. The Kalman filter (KF) is a popular method to online estimate unknown system state variables. This method conducts the state estimation in the form of feedback control, where system state variables are estimated using their measurements as feedback [27]. It can achieve the accuracy of minimum mean-square error with a small computational load. Comparing to other optimal filters such as the classical Wiener filter, KF works well for non-stationary signals and its estimation solution will not be distorted if the system state equation is constructed correctly [28]. However, KF can only estimate state variables in the time domain. On the other hand, FEM can solve state variables by decomposing the space domain into a finite number of elements. Therefore, the combination of KF and FEM provides a solution to estimate state variables not only in time but also in space. So far, there has been very limited research focusing on using Kalman filter for soft tissue deformation, and the combination of KF and FEM for modelling of soft tissue deformation is still a challenge task. Just recently, Yarahmadian et al studied a Kalman filter based on FEM for modelling of soft tissue deformation [29].
However, this work is in the preliminary stage, providing only simple simulation results on a cubic-shape model without any quantitative analysis. Further, since the system state equation is constructed from FEM modelling without considering system process noise, the resultant KF is not suitable for realworld situations where system process noise is inevitable. In addition, the use of the Newmark integration for time domain discretization cannot guarantee the stability for all circumstances.
It should also be noted that there has been very limited research on using the Kalman filter for soft tissue deformation, although the Kalman filter has been extensively studied in other areas such as signal processing, tracking and control, and vehicle navigation.
This paper presents a novel Kalman filter finite element method (KF-FEM) to predict soft tissue deformation in real time by combining KF with FEM. Rather than direct calculation of soft tissue deformation from FEM, which causes an expensive computational load, this method formulates FEM deformation of soft tissues as a filtering identification process to real-time estimate soft tissue behaviours based on local measurement data of displacement. It discretizes soft tissue deformation based on elastic theory in the space domain using FEM. Subsequently, the dynamic FEM model of soft tissue deformation is discretized in the time domain by the Wilson-θ implicit integration, and further transformed into a state space representation for filtering estimation. Based on above, a Kalman filter is developed to online estimate the dynamic deformation of soft tissues from local measurement of displacement. Interactive tool-tissue interaction with force feedback is also attained for surgery simulation. Simulations and experiments as well as comparison analysis with the traditional FEM have been conducted to comprehensively evaluate the performance of the proposed method. Different from the other methods improving the computational performance of FEM at the cost of losing the modelling accuracy, the proposed method significantly improves the computational performance but without sacrificing the modelling accuracy of FEM. This paper significantly improves and extends the authors' previous work reported in [29] with substantial further studies and development. It uses the Wilson-θ integration, which is unconditionally stable, to address the issue of conditional stability involved in the Newmark integration for time domain discretization. Further, the system state equation is constructed from FEM modelling by consideration of system process noise, making the resultant Kalman filter suitable for real-world situations. Consequently, additional and novel approaches are established for formulation of the proposed KF-FEM. In addition, this paper also provides comprehensive performance evaluation, including detailed quantitative analysis, in-depth simulation results and analysis on virtual human organs for surgery simulation, experimental verification and analysis, and haptic interaction with virtual human organs for surgical simulation. 53472 VOLUME 8, 2020

II. LINEAR ELASTICITY
Due to the complexity in structures and constituents, soft tissues exhibit complicated mechanical characteristics. Despite the complex material composition, the deformation of soft tissues can be accurately characterised by linear elasticity in a short time interval. The basic strain energy equation of linear elasticity can be defined as where X = [x, y, z] T represents the nodal position in the solution domain .
The relationship between the stress tensor σ and strain tensor ε is represented by Hooke's law of continuum media σ = Dε (2) where the matrix D is represented as and where E and v represent the Young's modulus and Poisson's ratio, respectively. Then the relationship between the strain and displacement can be written as ε = Bu (5) where the matrix B is represented as and N 1 , N 2 , N 3 , and N 4 represent the shape functions.

III. FINITE ELEMENT METHOD
Suppose that the soft tissue model is divided into a finite number of 4-node tetrahedral elements. The nodal displacement u = [u x , u y , u z ] T can be calculated as by interpolation Thus, the displacement gradient can be defined as From (1), the strain energy of an element can be written as The equilibrium equation of the element is given by where f e is the element force. Substituting (10) into (9) yields the following linear form where the stiffness matrix (K e ) for each finite element is (12) in which V e represents the element volume. Then, the global stiffness matrix and global force can be represented as

A. DYNAMIC EQUILIBRIUM SYSTEM
Assembling the equations of each finite element together and applying the mass and damping lead to the dynamic equilibrium system equation to represent the deformation behaviours of the soft tissue where the mass matrix M is a known function of material density; the stiffness matrix K is also a constant function of material constitutive law related to the material properties, i.e., Young's modulus E and Poisson's ratio v; and the stiffness and mass matrices are defined as The damping matrix C is frequency-dependent and determined by C = αM + βK with small weighting values of damping coefficient α and β [30]; f is the external force; and U represents the displacement,U the velocity andÜ the acceleration.

B. TEMPORAL DISCRETIZATION
The solution of the dynamic equilibrium system (15) in the time domain can be solved by the explicit or implicit integration. Despite its simplicity in implementation, the explicit integration is stable under small time steps only [3]. The implicit integration is superior to the explicit one for stiff equations. It can achieve the stability of dynamic simulation under large time steps, leading to the improved simulation efficiency [4]. The Wilson − θ is a popular implicit time integration scheme. Under the assumption that the acceleration varies linearly with time increment, this scheme uses an updated Lagrangian formulation to achieve a high numerical stability with a rapid convergence speed [31]. Given its advantages in both stability and convergence speed, the Wilson − θ scheme is adopted in this paper to solve the unknown displacement, velocity and acceleration in (14).
Using the Wilson − θ implicit time integration scheme, the acceleration and velocity at time t + θ t can be represented asÜ As shown in (16)-(18), to calculateÜ t+ t ,U t+ t and U t+ t , we need to know U t+θ t .
Applying the Wilson − θ implicit time integration scheme to (14) yields 6 Thus, we can readily have U t+θ t from (19). Since the Wilson − θ scheme is unconditionally stable when θ ≥1.37, in this paper the integration constant is set as θ = 1.4 [31].

IV. KALMAN FILTER FINITE ELEMENT METHOD A. STATE SPACE REPRESENTATION
In order to represent the state space model, the displacement, velocity and acceleration are considered as the unknown variables for soft tissue deformation. Rewrite (19) as (20) where I is the identity matrix, and A, B, Y and J are defined as It should be noted that the displacement calculated by (20) is for the time step θ t, rather than the actual time step t required in KF. To derive U t+ t for formulation of the state space equation, substituting (20) into (16) yields where A aa , A va , A xa and A fa are defined as Substituting (22) into (18) yields where A a , A v , A x and A f are defined as From (24) U t+ t can be calculated directly fromÜ t , U t , U t and f t .
Re-arranging (24) into state space representation gives where and where U 0 and P 0 denote the initial nodal displacement and its associated error covariance.

B. KALMAN FILTER
Based on the above state space representation by FEM modelling, the KF system state equation for dynamic soft tissue deformation is described as where the nodal displacement U t is the system state vector at time t, z t is the vector of control input at time t, F is the state transition matrix, G is the control input matrix, and w t is the process noise. The measurement equation of the dynamic soft tissue deformation system is defined as where y t is the measurement vector at time t, H is the measurement matrix, which is an identity matrix, and n t is the measurement noise. The process noise w t and measurement noise n t can be defined as where Q is the process noise covariance which is assumed to be a zero-mean with Gaussian distribution, and R is the measurement noise covariance which is assumed to be a zeromean with Gaussian distribution. It is also assumed that w t and n t are uncorrelated. The KF process involves two stages, which are the prediction and measurement update. The prediction stage is represented byŪ  whereŪ t+ t andP t+ t represent the priori estimated displacement and its associated error covariance. The measurement update stage is described as where K t+ t represents the Kalman gain. Since (33), (34) and (36) are independent of measurement, they can be calculated offline before the filtering algorithm starts. Fig. 1 shows the KF-FEM algorithm, which is made up of both offline (preprocessing) process and online (realtime computing) process.

V. PERFORMANCE ANALYSIS
A prototype system was implemented by the proposed KF-FEM for modelling of soft tissue deformation. This system was implemented on a PC platform with Intel R Core TM i7-8750 CPU @ 2.20 GHz, RAM 16.00 GB and 64-bit Windows 10. Simulations and experiments, together with comparison analysis with the traditional FEM, were conducted to comprehensively evaluate the performance of the proposed KF-FEM for soft tissue deformation. The displacements generated by the traditional FEM were taken as the true values for the comparison purpose. The estimation error was calculated by comparing the deformations estimated by the proposed KF-FEM with those by the traditional FEM. Interactive tooltissue interaction with haptic feedback was also examined for the proposed KF-FEM.

A. CUBE MODEL
Simulation analysis was carried out on a cube model in the size of 6 × 6×6 mm 3 (see Fig. 2) by both FEM and KF-FEM. This cube model was uniformly meshed with 413 nodes and 1723 tetrahedral elements. The mechanical properties of the human liver [5], [32], i.e., mass density = 1060kg/m 3 , Young's modular = 3500Pa and Poisson's ratio = 0.49, were assigned to this model. A mechanical test was conducted on the cube model, where the bottom face was fixed, and a linear tensile force varied from 0 N to 0.001 N in the z direction was applied to the top face. The damping coefficients α and β were set as 0.05 and 0.0. The initial state value was set to 0. The process and measurement noise covariances were set to 0.01 and 0.001. The total simulation time was 5 s, and the time step ( t) was 0.02 s. To conduct an in-depth analysis on the performance of KF-FEM, two nodes in the cube model (except for the bottom face) were selected to observe the deformation displacements. One is the center point of the top face (observation node A in Fig. 2), where the measurement data of displacement, as shown in Fig. 3, were acquired by adding a random Gaussian white noise in the true values (i.e., the displacements by the traditional FEM). The other is a random node without the applied force (i.e., observation node B in Fig. 2). Fig. 4 shows the deformation results by both FEM and KF-FEM for the tensile test of the cube model. Both results are almost identical without any visually identified difference. To further investigate the difference between the KF-FEM and FEM results, Fig. 5 shows the displacements at the two observation nodes. It can be seen that the displacements at both observation nodes by the proposed KF-FEM are very close to those by FEM. As shown in Table 1, the mean error, RMSE (root mean square error) and standard deviation of KF-FEM are 5.90 × 10 −5 mm, 0.0015756 mm and 0.0015776 mm for observation node A, and 1.55 × 10 −5 mm, 0.0015709 and 0.0015739 mm for observation node B. Further, as shown in Table 2, the computational performance of KF-FEM is about 16 times faster than that of FEM, where the total computational time at each iteration is 0.05676 s for FEM, while 0.00348 s for KF-FEM.
It is evident from the above results that the proposed KF-FEM clearly filters the measurement noises shown in Fig. 3, leading to the similar accuracy as the traditional FEM. Further, the proposed KF-FEM also overcomes the FEM limitation of expensive computations. The FEM refresh rate for both tensile and compression tests of the cube model   is about 17 Hz. This means that FEM is unable to satisfy the minimum requirement of real-time visual feedback, i.e., 30 Hz refresh rate. In comparison, the KF-FEM refresh rate is about 280 Hz, which is much higher than the minimum requirement for real-time visual feedback. Thus, the proposed KF-FEM achieves the real-time performance but without decreasing the physical realism of FEM for soft tissue deformation.

B. HUMAN LIVER MODEL
Simulation analysis was also carried out on a virtual human liver model (see Fig. 6) by both FEM and KF-FEM. This liver model was uniformly meshed with elements with 1083 nodes and 4941 tetrahedral elements. The material properties of the human liver as above were also assigned to the virtual liver model. A mechanical test was conducted on the liver model, where the bottom surface of the liver model (indicated by the black nodes) was assumed fixed (see Fig. 6(a)), and a linear compression force varied from 0 N to 0.003 N was applied in the z direction to a local area (indicated by the red nodes) on the top surface of the liver model (see Fig. 6(b)). The other simulation parameters were the same as the case of the cube model.
Similar to the case of the cube model, two observation nodes were selected in the liver model (except for the bottom surface) to evaluate the deformation accuracy. As shown in Fig. 6(b), observation node A was the node where the measurement data of displacement were acquired, and observation node B was selected from the nodes without the applied compression force. Fig. 7 shows the noisy measurement data acquired at observation node A by adding a random Gaussian white noise in the true values for the compression test. Fig. 8 illustrates the deformations of the liver model by both FEM and KF-FEM for the compression test. It is evident that the deformation results of both FEM and KF-FEM are almost identical without any visually identified difference. To look into the difference between the FEM and KF-FEM deformation results, Fig. 9 shows the displacements of both FEM and KF-FEM at the two observation nodes for the compression test of the liver model. It can be seen that the proposed KF-FEM clearly filters the measurement noise shown in Fig. 7 and its displacements at both observation nodes are also very close to those of the traditional FEM. As shown in Table 3, the mean error, RMSE and standard deviation by the proposed KF-FEM are 7.85 × 10 −5 mm, 0.0017267 mm and 0.0017284 mm for observation node A, and 7.52 × 10 −5 mm, 0.0017251 mm and 0.0017269 mm for observation node B.
The computational performance of KF-FEM for the compression test of the liver model is about 25 times faster than that of FEM. As shown in Table 4, the computational time at each iteration is 0.48 s for FEM, while 0.01937 s for KF-FEM. The visual refresh rate of FEM is about 2.08 Hz, which is far to meet the minimum requirement (30 Hz refresh rate) for real-time visual feedback. In contrast, the visual refresh rate of KF-FEM is about 52 Hz, which is about 26 times larger than that of FEM and is sufficient to achieve real-time visual feedback. Thus, it is evident that the proposed KF-FEM overcomes the FEM limitation of expensive computations and satisfies the requirement of real-time performance for soft tissue simulation.
The above simulation results demonstrate that the proposed KF-FEM not only overcomes the FEM limitation of   expensive computations, but also inherits the FEM advantage of high modelling accuracy for soft tissue deformation.

C. EXPERIMENTAL VERIFICATION
Experimental trials were also conducted to evaluate the performance of the proposed KF-FEM for dynamic modelling of soft tissue deformation based on real measurement data of displacement. As shown in Fig. 10, a compression test was conducted on a phantom soft tissue sample using the robotic indentation system [33] to acquire the measurement data of displacement. The phantom tissue was made up of silicone rubbers (Young's modular = 50 kPa and Poisson's ratio = 0.49) with the similar characteristics of human soft tissues [34]. It was in cubic shape with the size of 10 × 10 × 5 mm. Its density was 1072 kg m 3 . Its bottom face was fixed on the blue solid platform. A frictionless rigid conic-tip probe with the flat end of 1 mm diameter was VOLUME 8, 2020 FIGURE 11. The cube model is the same as the phantom tissue sample in terms of geometry, material properties and loading conditions with the fixed boundary condition on the bottom face (indicated by the black and gray nodes). The force was applied on the top face (indicated by the red points), where the measurement data of displacement were acquired at observation node A and one random node without measurement acquired from (observation node B) was also selected to observe the deformation displacements by both FEM and KF-FEM. used to compress the phantom tissue in the vertical direction. A linear compression force from 0 N to 0.14 N was applied to the center of the phantom tissue's top face in the time period of 3.5 s.
For comparison analysis, both FEM and KF-FEM modelling under the same circumstances as the experimental compression test were performed, where the virtual tissue model (see Fig.11), which was uniformly discretized with 826 nodes and 3629 tetrahedral elements, was the same as the phantom tissue sample in terms of geometry, material properties and loading conditions. The measurement data of displacement at the contact node where the compression force was applied were collected from the experimental compression test on the phantom tissue sample. Based on the acquired measurement data, the deformations of the phantom tissue were estimated by the proposed KF-FEM, and further compared with the deformations by the traditional FEM. Fig. 12 shows the deformation results of both FEM and KF-FEM are in an excellent agreement. To further investigate the deformation difference between FEM and KF-FEM for the experimental test, in addition to the contact node (observation node A) where the measurement data of displacement were acquired, one more random node (observation node B) in the phantom tissue model (except for the bottom face) was selected to observe the deformations by both FEM and KF-FEM. As shown in Fig. 13, the proposed KF-FEM clearly filters the noisy measurement data, and its estimated displacements are well in line with those by the traditional FEM. As shown in Table 5, the mean error, RMSE and standard deviation of the proposed KF-FEM are 0.0011318 mm, 0.033626 mm and 0.033987 mm, for observation node A, and 0.0001455 mm, 0.001158 mm and 0.001161 mm for observation node B. Fig. 14 shows the computational performances of both FEM and KF-FEM for the experimental compression test. It can be seen that the traditional FEM achieves the refresh rate of 30 Hz at about 1200 elements, while the proposed KF-FEM at about 5300 elements. This proves the presented KF-FEM is much superior to the traditional FEM in terms of computational time. It should also be noted from the above simulation and experimental analyses that the proposed KF-FEM estimates the deformation on the entire model from the only measurement data of displacement at the contact node where the external force is applied. In other words, the proposed KF-FEM estimates the global deformation of the entire model from local measurement of displacement at the contact node where the external force is applied.

D. FORCE FEEDBACK
Tool-tissue interaction with force feedback was also attained by integration of the presented KF-FEM with a PHANToM  haptic instrument from Geomagic. The measurement data of displacement were acquired from the penetration depth of PHANToM's stylus tip which is in contact with the virtual tissue model. The elastic force was calculated by the proposed KF-FEM and transferred to the haptic instrument for force feedback. Fig. 15 shows the deformations of the virtual human liver model in Section V-B via the virtual haptic needle.
To achieve realistic force feedback by the PHANToM haptic device, the force refresh rate must be at least 1,000 Hz. As shown in Fig. 16, the proposed KF-FEM achieves the force refresh rate of 1,000 Hz at around 1100 tetrahedron elements, while the traditional FEM at around 100 tetrahedron elements. This demonstrates that the proposed KF-FEM is much superior to the traditional FEM in terms of haptic performance. In the case that the required haptic refresh rate cannot be met, the force information at current iteration is extrapolated from the previous iteration for realistic force feedback [4], [35].

VI. CONCLUSION
This paper proposes a new KF-FEM for dynamic modelling and analysis of soft tissue deformation. It carries out soft tissue deformation with a filtering identification process to online estimate mechanical behaviours of soft tissues based on local measurement of displacement. The contribution of this paper is that a Kalman filter and its associated system equations are established for online estimation of soft tissue deformation. Soft tissue deformation is discretized in space by FEM and in time by the Wilson − θ implicit integration to formulate the system state equation for filtering estimation. Based on the system state equation, a Kalman filter is developed to provide the statistical estimation of soft tissue deformation from local measurement of displacement. Tooltissue interaction with force feedback is also attained for surgery simulation. Simulations and experiments as well as comparison analysis demonstrate that the proposed KF-FEM significantly improves the FEM computing performance but without sacrificing the FEM modelling accuracy for soft tissue deformation. The proposed KF-FEM also improves the simulation efficiency by taking a large time step for dynamic modelling of soft tissue deformation.
Future research work will focus on improvement of the proposed method. Currently, the proposed method considers linear elastic behaviours only. In the future, based on nonlinear elasticity, the linear system state equation will be extended to formulate a nonlinear state equation to characterize nonlinear mechanical behaviours of soft tissues. Consequently, a nonlinear filtering algorithm such as the extended Kalman filter [36] or unscented Kalman filter [37] will be developed based on the nonlinear system state equation to real-time estimate nonlinear deformation behaviours of soft tissues.