Soft Tissue Deformation Modeling in the Procedure of Needle Insertion: A Kriging-based Method

The simulation and planning system (SPS) requires accurate and real-time feedback regarding the deformation of soft tissues during the needle insertion procedure. Traditional mechanical-based models such as the ﬁnite element method (FEM) are widely used to compute the deformations of soft tissue. However, it is diﬃcult for the FEM or other methods to ﬁnd a balance between an acceptable image ﬁdelity and real-time deformation feedback due to their complex material properties, geometries and interaction mechanisms. In this paper, a Kriging-based method is applied to model the soft tissue deformation to strike a balance between the accuracy and eﬃciency of deformation feedback. Four combinations of regression and correlation functions are compared regarding their ability to predict the maximum deformations of ten characteristic markers at a ﬁxed insertion depth. The results suggest that a ﬁrst order regression function with Gaussian correlation functions can best ﬁt the results of the ground truth. The functional response of the Kriging-based method is utilized to model the dynamic deformations of markers at a series of needle insertion depths. The feasibility of the method is veriﬁed by investigating the adaptation to step variations. Compared with the ground truth of the ﬁnite element (FE) results, the maximum residual is less than 0 . 92 mm in the Y direction and 0 . 31 mm in the X direction. The results suggest that the Kriging metamodel provides real-time deformation feedback for a target and an obstacle to a SPS.

In needle-based minimally invasive surgeries, needle insertion errors always lead to a  [27]. It is a simplified computer-based simulation model that has an input/output function based on the surface response [28]. The advantages of 58 the Kriging metamodel are twofold: First, the Kriging metamodel can reduce the 59 large requirement for computing resources and offer real-time simulation, especially 60 being able to overcome the challenges of computational techniques in complex simu-61 lation models [29]. Hence, depending on the accurate simulation dataset generated 62 by other offline needle-tissue interaction models, the online Kriging model can be 63 established rapidly without reducing the accuracy of the original model by much. 64 Second, the computer-based simulation program is deterministic, that is, the same 65 input corresponds to the same output, which cannot fully reflect the uncertainty of 66 mechatronic systems [30,31]. In robot-assisted needle insertion systems, complex 67 factors such as the imaging equipment, physician's skill and patient's condition lead 68 to uncertainty regarding tissue deformation. The Kriging metamodel with the ran-69 dom sampling method can reflect the uncertainties of robot-assisted needle-tissue 70 interactions. 71 In our work, the input dataset is generated by Latin hypercube sampling (LHS), 72 the parameters of which include the material properties of the tissues and needle, the 73 geometrical properties of the needle and the solver parameters. Ten characteristic 74 markers are selected to represent the targets and obstacles inside the tissue body.

75
The output dataset is generated by running the needle-tissue coupling FE model 76 presented in [32], which is used as the ground truth of the Kriging-based model.

77
The validity and feasibility of the proposed Kriging-based model are analyzed, and 78 the results suggest that the combination of the Kriging metamodel and the high-79 precision finite element model provides real-time deformation feedback for a target 80 and an obstacle to the SPS.

81
The rest of this paper is organized as follows. Section II introduces the data   Kriging-based prediction model, ten characteristic markers in the FE model are  There were 11 variables chosen and listed in Table 1, wherein x1, x2, and x11 119 are the material parameters of soft tissues, x3 -x7 are the needle parameters, x8 120 -x9 are the FEM solver parameters and x10 is the friction coefficient between the 121 needle and soft tissue.

122
The input variables are denoted as X = [x 1 , · · · , x 11 ], which are employed to construct the design matrix of the Kriging-based model. Running the finite element model 20 times, the displacements at the 10 observed locations are collected as the output responses of the Kriging model, denoted by Y = u 1 xy , · · · , u 10 xy T = [y 1 , y 2 , · · · , y 10 ] T , y i ∈ R 20 . The soft tissue deformation occurs in both the x and y directions in the two-dimensional case; hence, the output response u xy is the resultant displacement in both directions and is written as  where u x and u y are the displacements in the x and y directions, respectively.

123
Since the resultant displacement increases as the needle insertion depth into the 124 tissue increases, the maximum value of u xy is chosen as the output response of the 125 Kriging-based prediction model.

126
Kriging model construction

127
The Kriging prediction model is an interpolation of known observation locations, of which the mean square error equals zero. The Kriging model is usually expressed as a combination of a polynomial and its deviation [33], written aŝ where f (x) is the basis polynomial function, usually a regression function. r(x) is the column vector of correction matrix where F is the design matrix, written as m is the number of groups of FE model, and in this application, m = 20. p is the number of basis functions, which is determined by the type of basis function. In our experiments, p = 12 for first order regression functions. The variable x i is obtained from the FE model by the Latin hypercube sampling method, and x i is normalized to the interval [0, 1], as iñ wherex i is the normalized data of x i and x L and x U are the minimum and maximum 128 values of the variable x i , respectively.

129
The mean square error of the Kriging prediction model is formulated as (5) [34]

132
The design matrix F m×p in the Kriging prediction model is where the basis function f (x) has several forms. Usually, constant and linear forms 133 of the regression function are used to construct the basis function [35].

134
• Constant form, p = 1, of which the design matrix F is a column vector.
• Linear form, p = n + 1, of which the basis function is written as (8).
where the component form of the estimate point x is expanded as x = [w 1 , w 2 , · · · , w n ] ∈ R n . In our experiment, both zero order (constant form, p=1) and first order (linear 136 form, p=12) functions are adopted and compared.

Correlation function 138
The Kriging model assumes that the correlation of the output is determined by the distance between the input variables. The correlation function of the input variables is written as the product of n one-dimensional correlation equations, as shown in (9): where d j = w j − x j and θ is the correlation coefficient. In the Kriging model, the most widely used correlation model is that shown in (10): where · denotes the Euclidean distance of d. When p = 1 and p = 2, the exponential and Gaussian correlation functions are as shown in (11) and (12).
The correlation function decreases with the increase in the Euclidean distance d j , and a larger correlation coefficient θ leads to a rapid decline of the correlation function. Substituting the correlation function into (5), the mean square error is a function of the covariance σ 2 and the correlation coefficient θ. The optimal solution to the correlation coefficient θ * is converted into an unconstrained global optimization problem [36], as shown in (13): where |R| is the determinant value of the correlation matrix R. The input and output data are prepared by running the finite element program 20 times, as in [32]. To establish a Kriging model, 19 sets of observation points ([X in (1, : ), · · · , X in (i − 1, :), X in (i + 1, :), · · · , X in (20, :)]) are chosen as the input design sites, and the remaining set X in (i, :) is used as the test set. The zero order and first order regression functions and the Gaussian and exponential correlation functions are employed to construct the Kriging model for static soft tissue deformation. In this section, the maximum displacements of the marked points are regarded as the static soft tissue deformations. Therefore, there are four kinds of Kriging model, namely, zero order Gaussian, zero order exponential, first order Gaussian and first order exponential, which are shown in Fig. 3. In the figure, the experimental values  The test set of each run is randomly selected. The upper-left test set is X in (5, :), the upper-right test set is X in (6, :), the lower-left test set is X in (7, :), and the lower-right test set is X in (8, :). From Fig. 3, it is also suggested that the closer the tissue is to the needle body, the greater the deformation. To analyze the impacts of the regression function and correlation function on the model, the residual e of the simulation is defined as the difference between the finite element simulation value y and the Kriging prediction valueŷ, that is, e =ŷ − y. The maximum residual values of the zero order Gaussian, zero order exponential, first order Gaussian and first order exponential Kriging models are 7.5108 mm, 7.3327 mm, 3.2752 mm and 3.2994 mm, respectively. Furthermore, the mean residual of each output response is defined as in (14): where m is the number of finite element program samples, i.e., m = 20. i is the position number of the marked positions, i = 1, · · · , 10. The mean residuals of the 10 observed locations are shown in Fig. 4. In Fig. 4, the mean residuals of the zero order Gaussian and exponential types are 1.5650 mm and 1.3992 mm, respectively. The mean maximum residuals of the first order Gaussian and exponential are 0.8941 mm and 0.9177 mm, respectively. From the prediction results, the first order regression function is better than the zero order regression function, and the Gaussian correlation function is slightly better than the exponential function from 1, · · · , 11). When the output displacement increases with the variables, the SI is 157 greater than zero. In the opposite case, the SI is less than zero. According to the 158 range of the SI, the 11 input variables can be classified into three types, i.e., leading 159 parameters, nonleading parameters and disturbed parameters, as shown in Fig. 5,   160 where the abscissa represents 11 variables and the ordinate represents the SI. Since where x = x i + y j + z k is the position vector at the initial time, X = X i + Y j + Z k 199 is the position vector at time t, and u = u x i + u y j + u z k is the displacement of

208
The input variables in the FE simulation experiment are x = [x 1 , x 2 , · · · , x n ] T , x i ∈ R n , i = 1, 2, · · · , m, and the output response is a function of time t. The output response where r i denotes the time step. The Kriging model of the functional response is written as (17) on the basis of (2): where y(x, t) is the output response of the input variable x at time t. F(x, t) =

209
[f 1 (x, t), · · · , f p (x, t)] T is a set of known polynomial basis functions, where usually 210 f 1 (x, t) = 1. β is an unknown basis function coefficient. z(x, t) is a zero-mean 211 Gaussian random function, the covariance function of which is shown in (18): Assume that the correlation function R(x 1 − x 2 , t 1 − t 2 ) is the product of n onedimensional correlation equations, as shown in (19): where R i (x i1 − x i2 ) is the correlation function of the i th set of variables and R T (t 1 − 213 t 2 ) is the correlation function at time t.
where ⊗ is the Kronecker product operator and 1 ri is the column vector with length r i = 1. When the time steps are the same, X is an mr × n matrix.
is the corresponding functional space.
Combined with equations (2) and (17), the Kriging model of the functional response is written as shown in (21): where F is written as in which R X,t is an N × N correlation matrix, the elements of which are R(X i − 215 X j , t * i − t * j ). The correlation coefficient θ is optimized with (13). The optimization 216 of the correlation coefficient requires a large number of solutions of R −1 X,t and |R X,t | 217 so that the computational cost is extremely high.

218
In this paper, we consider the functional output on the regular grid, i.e., t 1 =· · ·= t n = t , r 1 = · · · = r n = r. The correlation matrix R X,t of the functional response is written as (23): where R X is an m×m-dimensional correlation matrix whose elements are R(x 1 −x 2 ) and R t is an r×r-dimensional correlation matrix whose elements are R T (t i −t j ). The Gaussian correlation function is used to construct the correlation matrix, and the inverse of the correlation matrix is simplified to R −1 X,t = R −1 X ⊗ R t −1 . Substituting this matrix into (21), we can write the Kriging model of the functional response as (24):ŷ The computational complexity of the inversion matrix is reduced from O(n 3 m 3 ) to 219 O(n 3 + m 3 ) with the Kronecker product operator.
(17) takes account of the coupling effect between the input variable x and time t. The basis function F(x, t) is identified by estimating the average value of the functional response. In the case of r 1 =· · · =r n =r, we defineē ·j = 1 m m i=1 (y ij −ȳ i· ) andȳ i· = 1 r r l=1 y il .ē = [ē ·1 , · · · ,ē ·r ] T andȳ = [ȳ 1· , · · · ,ȳ m· ] are fitted as shown in (25).ē The unknown polynomial basis functions k(x) and g(x) are constructed us-221 ing the regression function. Furthermore, the basis function in (24) is writ- where r is the total number of puncturing steps, i.e., r = 31, and i is the position 226 number of the output response, where i = 1, · · · , 10.   Table 2. Overall, the Kriging model smooths the original finite element *·max maximum residual. *·mean mean residual. *·re relative mean residual.   nine observation positions are listed in Table 3. From Table 3, the maximum resid-277 ual in the X direction occurs at N 1 and N 7; its value is 0.31 mm. The maximum 278 residual in the Y direction appears at N2; its value is 0.92 mm. The relative mean 279 residual in the Y direction is less than 8% and has a relatively high accuracy. For 280 the real insertion procedure, the displacements in the X direction are very small, 281 sometimes approaching zero; hence, the relative mean residual is larger than 20%.

282
However, the absolute mean residual of each location is less than 0.31 mm, which 283 can meet the estimation accuracy.  The performance of the proposed Kriging-based model depends on the accuracy 303 of the dataset generated by the mechanical-based simulation greatly. Therefore, 304 future work will focus on combinations with a more accurate simulation model, 305 considering complex tissue behavior ranging from hyperelasticity to viscoelasticity, 306 and the needle deflection will be modeled at the same time. The sensitivity of 307 additional parameters will be analyzed to improve the accuracy and efficiency of 308 the Kriging-based model. The datasets analysed during the current study are available from the corresponding 312 author on reasonable request.