# Soft Tissue Deformation Modeling in the Procedure of Needle Insertion: A Kriging-Based Method – Chinese Journal of Mechanical Engineering

Sep 5, 2022

### Kriging Model of the Functional Response

According to the continuum mechanics, the motion equation of a point inside the tissue (mathbf{O}) in the Euclidean space is expressed as Eq. (15):

$${varvec{X}}=chi left({varvec{x}},tright)mapsto u={varvec{X}}-{varvec{x}}=chi left({varvec{x}},tright)-{varvec{x}},$$

(15)

where ({varvec{x}}=xoverrightarrow{i}+yoverrightarrow{j}+zoverrightarrow{k}) is the position vector at the initial time, ({varvec{X}}=Xoverrightarrow{i}+Yoverrightarrow{j}+Zoverrightarrow{k}) is the position vector at time t, and ({varvec{u}}={u}_{x}overrightarrow{i}+{u}_{y}overrightarrow{j}+{u}_{z}overrightarrow{k}) is the displacement of point (mathbf{O}) at time t. The displacement of tissue nodes is related to the time and space coordination, and it is a functional response with time and space indices. When the Kriging metamodel is used to solve the functional response problem, the calculation of the correlation coefficient of the maximum likelihood estimation includes the inverse matrix and the determinant of the correlation matrix, which leads to the instability and time complexity of the computation [35]. The Kronecker product is employed to construct the correlation matrix and to reduce the computational burden of the Kriging method [36].

The input variables in the FE simulation experiment are ({varvec{x}}={[{x}_{1},{x}_{2},cdots ,{x}_{n}]}^{mathrm{T}}), ({{varvec{x}}}_{i}in {mathbf{R}}^{n},i=mathrm{1,2},cdots ,m), and the output response is a function of time (t). The output response ({{varvec{y}}}_{i}={[{y}_{{i}_{1}},{y}_{{i}_{2}},cdots ,{y}_{{i}_{{r}_{i}}}]}^{mathrm{T}}) is recorded at each time ({{varvec{t}}}_{i}=[{t}_{{i}_{1}},{t}_{{i}_{2}},cdots ,{t}_{{i}_{{r}_{i}}}{]}^{mathrm{T}}), where ({r}_{i}) denotes the time step. The Kriging model of the functional response is written as Eq. (17) on the basis of Eq. (2):

$$y({varvec{x}},t)=mathcal{F}({varvec{x}},t)beta +z({varvec{x}},t),$$

(16)

where (y({varvec{x}},t)) is the output response of the input variable ({varvec{x}}) at time t. (mathcal{F}({varvec{x}},t)=[{f}_{1}({varvec{x}},t),cdots ,{f}_{p}({varvec{x}},t){]}^{mathrm{T}}) is a set of known polynomial basis functions, where usually ({f}_{1}({varvec{x}},t)=1). (beta) is an unknown basis function coefficient. (z({varvec{x}},t)) is a zero-mean Gaussian random function, the covariance function of which is shown in Eq. (17):

begin{aligned} Vleft( {{varvec{x}}_{1} ,{varvec{x}}_{2} } right) =& covleft( {zleft( {{varvec{x}}_{1} ,t_{1} } right),,zleft( {{varvec{x}}_{2} ,t_{2} } right)} right. \ = & sigma^{2} Rleft( {{varvec{x}}_{1} – {varvec{x}}_{2} ,t_{1} – t_{2} } right). \ end{aligned}

(17)

Assume that the correlation function (R({{varvec{x}}}_{1}-{{varvec{x}}}_{2},{t}_{1}-{t}_{2})) is the product of (n) one-dimensional correlation equations, as shown in Eq. (18):

$$R({{varvec{x}}}_{1}-{{varvec{x}}}_{2},{t}_{1}-{t}_{2})=left[prod_{j=1}^{n}{R}_{i}({x}_{i1}-{x}_{i2})right]{R}_{T}({t}_{1}-{t}_{2}),$$

(18)

where ({R}_{i}({x}_{i1}-{x}_{i2})) is the correlation function of the ith set of variables and ({R}_{T}({t}_{1}-{t}_{2})) is the correlation function at time (t).

The functional response is assumed to be the ((Ntimes 1))-dimensional vector ({varvec{Y}}=[{{varvec{y}}}_{1}^{mathrm{T}},{{varvec{y}}}_{2}^{mathrm{T}},cdots ,{{varvec{y}}}_{n}^{mathrm{T}}{]}^{mathrm{T}}), and (N=sum_{i=1}^{m}{r}_{i}) . When the time step is (r), then (N=mr). The corresponding ((Ntimes n))-dimensional design matrix ({varvec{X}}=[{{varvec{X}}}_{1},{{varvec{X}}}_{2},cdots ,{{varvec{X}}}_{N}]) is represented with the Kronecker notation as shown in Eq. (19):

$$begin{array}{c}begin{array}{cc}{{varvec{X}}}_{Ntimes n}& ={ left[{1}_{{r}_{1}}^{mathrm{T}}otimes {{varvec{x}}}_{1},cdots ,{1}_{{r}_{n}}^{mathrm{T}}otimes {{varvec{x}}}_{n}right]}^{mathrm{T}}\ & =left[begin{array}{cccc}{x}_{11}& {x}_{12}& cdots & {x}_{1n}\ vdots & vdots & vdots & vdots \ {x}_{11}& {x}_{12}& cdots & {x}_{1n}\ {x}_{21}& {x}_{22}& cdots & {x}_{2n}\ vdots & vdots & vdots & vdots \ {x}_{m1}& {x}_{m2}& cdots & {x}_{mn}end{array}right]begin{array}{c}mapsto {text{R}}{text{o}}{text{w}} , {1}\ vdots \ mapsto {text{R}}{text{o}}{text{w}} , {r}_{i}\ mapsto {text{R}}{text{o}}{text{w}} , {r}_{i}text{+}{1}\ vdots \ mapsto {text{R}}{text{o}}{text{w}} , {N}end{array},end{array}end{array}$$

(19)

where (otimes) is the Kronecker product operator and (varvec{1}_{{r}_{i}}) is the column vector with length ({r}_{i}={1}). When the time steps are the same, ({varvec{X}}) is an (mrtimes n) matrix. ({varvec{T}}=[{{varvec{t}}}_{1}^{mathrm{T}},cdots ,{{varvec{t}}}_{n}^{mathrm{T}}{]}^{mathrm{T}}={[{{varvec{t}}}_{1}^{*},{{varvec{t}}}_{2}^{*},cdots ,{{varvec{t}}}_{N}^{*}]}^{mathrm{T}}in {mathbf{R}}^{Ntimes 1}) is the corresponding functional space. Combined with Eq. (2) and Eq. (16), the Kriging model of the functional response is written as shown in Eq. (20):

$$hat{y}({varvec{x}},{varvec{t}})=F({varvec{x}},{varvec{t}})hat{{varvec{beta}}}+{R}^{mathrm{T}}({varvec{x}},{varvec{t}}){R}_{{varvec{X}},{varvec{t}}}^{-1}({varvec{Y}}-{varvec{F}}hat{{varvec{beta}}}),$$

(20)

where F is written as

begin{aligned} F & = {[}F{(}X_{{1}} {,};t_{1}^{ * } {),} cdots {,};F{(}X_{N} {,};t_{N}^{ * } {)]}^{{text{T}}} {,} \ hat{user2{beta }} & = {(}{varvec{F}}^{{text{T}}} {varvec{R}}_{{{X}{,},{t}}}^{ – 1} {varvec{F}}{)}^{ – 1} {varvec{F}}^{{text{T}}} {varvec{R}}_{{X{,},t}}^{ – 1} {mathbf{Y}}{,} \ R{(}{varvec{x}}{,};{varvec{t}}{)} & { = [}R{(}{varvec{x}} – {varvec{X}}_{{1}} {,};{varvec{t}} – {varvec{t}}_{1}^{ * } {),} cdots {,};R{(}{varvec{x}} – {varvec{X}}_{N} {,};{varvec{t}} – {varvec{t}}_{N}^{ * } {)]}^{{text{T}}} {,} \ end{aligned}

(21)

where ({varvec{R}}_{{{X}{,},{t}}}^{{}}) is an (Ntimes N) correlation matrix, the elements of which are (R({{varvec{X}}}_{i}-{{varvec{X}}}_{j},{{varvec{t}}}_{i}^{*}-{{varvec{t}}}_{j}^{*})). The correlation coefficient (theta) is optimized with Eq. (13). The optimization of the correlation coefficient requires a large number of solutions of ({varvec{R}}_{{{X}{,},{t}}}^{ – 1}) and (left| {{varvec{R}}_{{{X}{,},{t}}}^{{}} } right|) so that the computational cost is extremely high.

In this paper, we consider the functional output on the regular grid, i.e., ({{varvec{t}}}_{1}=cdots ={{varvec{t}}}_{n}={varvec{t}}), ({{varvec{r}}}_{1}=cdots ={{varvec{r}}}_{n}={varvec{r}}). The correlation matrix ({varvec{R}}_{{{mathbf{X}}{,},{mathbf{t}}}}^{{}}) of the functional response is written as Eq. (22):

$${{varvec{R}}}_{{varvec{X}},{varvec{t}}}={{varvec{R}}}_{{varvec{X}}}otimes {{varvec{R}}}_{{varvec{t}}},$$

(22)

where ({{varvec{R}}}_{{varvec{X}}}) is an ((mtimes m))– dimensional correlation matrix whose elements are ({varvec{R}}({{varvec{x}}}_{1}-{{varvec{x}}}_{2})) and ({{varvec{R}}}_{{varvec{t}}}) is an ((rtimes r))-dimensional correlation matrix whose elements are ({{varvec{R}}}_{{varvec{T}}}({{varvec{t}}}_{i}-{{varvec{t}}}_{j})). The Gaussian correlation function is used to construct the correlation matrix, and the inverse of the correlation matrix is simplified to:

({{varvec{R}}}_{{varvec{X}},{varvec{t}}}^{-1}={{varvec{R}}}_{{varvec{X}}}^{-1}otimes {{{varvec{R}}}_{{varvec{t}}}}^{-1}). Substituting this matrix into Eq. (20), we can write the Kriging model of the functional response as Eq. (23):

$$hat{y}({varvec{x}},{varvec{t}})={varvec{F}}({varvec{x}},{varvec{t}})hat{{varvec{beta}}}+{{varvec{R}}}^{mathrm{T}}({varvec{x}},{varvec{t}}){{varvec{R}}}_{{varvec{X}}}^{-1}otimes {{{varvec{R}}}_{{varvec{t}}}}^{-1}({varvec{Y}}-{varvec{F}}hat{{varvec{beta}}}).$$

(23)

The computational complexity of the inversion matrix is reduced from (O({n}^{3}{m}^{3})) to (O({n}^{3}+{m}^{3})) with the Kronecker product operator.

Eq. (16) takes account of the coupling effect between the input variable ({varvec{x}}) and time (t). The basis function (F({varvec{x}},t)) is identified by estimating the average value of the functional response.

In the case of ({r}_{1}=cdots ={r}_{n}=r), we define ({overline{{varvec{e}}}}_{cdot j}=frac{1}{m}sum_{i=1}^{m}({y}_{ij}-overline{{y}_{icdot }})) and ({overline{y}}_{icdot }=frac{1}{r}sum_{l=1}^{r}{y}_{il}). (overline{{varvec{e}}}=[{overline{e}}_{cdot 1},cdots ,{overline{e}}_{cdot r}{]}^{mathrm{T}}) and (overline{{varvec{y}}}=[{overline{y}}_{1cdot },cdots ,{overline{y}}_{mcdot }]) are fitted as shown in Eq. (24):

begin{aligned} overline{user2{e}}left( t right) & = k^{{text{T}}} left( t right)beta_{t} + zleft( t right), \ overline{user2{y}}left( {varvec{x}} right) & = g^{{text{T}}} left( {varvec{x}} right)beta_{x} + zleft( {varvec{x}} right). \ end{aligned}

(24)

The unknown polynomial basis functions (k(x)) and (g(x)) are constructed using the regression function. Furthermore, the basis function in Eq. (23) is written as ({varvec{F}}({varvec{x}},t)=[1,{k}^{mathrm{T}}(x),{g}^{mathrm{T}}(x){]}^{mathrm{T}}) and the matrix ({varvec{F}}) is expressed as ({varvec{F}}=[F({X}_{1},{t}_{1}^{*}),cdots ,F({X}_{N},{t}_{N}^{*}){]}^{mathrm{T}}).

### Deformation Modeling of Needle Insertion at a Series of Depth

One hundred sets of input variables are generated with LHS to construct the functional response Kriging model. After the normalized inverse operation, the finite element program runs 100 times to generate the deformations in the (X) and (Y) directions at 10 marked points for 31 time steps, as shown in Figure 2. For needle insertion at a constant speed, the insertion time is replaced with the depth (d). Here, the 31 time steps are expressed as ({varvec{d}}=[1.0text{ mm}, 2.0text{ mm}, 5.0text{ mm},cdots , 89text{ mm}]), and the increment of the insertion depth is (text{3 mm}), except in the first step.

Forty sets of data are randomly selected from 100 sets of records and used to construct the functional Kriging prediction model of soft tissue deformation, and an independent set of data is used to test the model. According to the former analysis, the polynomial basis functions are first order regression functions, and the correlation functions are Gaussian random functions. To evaluate the real-time performance of the Kriging module, the mean residual ({overline{widetilde{{varvec{e}}}}}_{i}) and the relative mean residual ({overline{widetilde{varepsilon }}}_{i}) at each time step (insertion depth) are defined for each observation location, written as Eq. (25):

begin{aligned} overline{{tilde{user2{e}}}}_{{i}} & = frac{1}{{r}}mathop sum limits_{j = 1}^{{r}} |{varvec{e}}_{{ji}} left| { = frac{1}{{r}}mathop sum limits_{j = 1}^{{r}} |hat{y}_{{ji}} – y_{{ji}} } right|, \ overline{{tilde{varepsilon }}}_{{i}} & = frac{1}{{r}}mathop sum limits_{j = 1}^{{r}} frac{{left| {hat{y}_{{ji}} – y_{{ji}} } right|}}{{left| {y_{{ji}} } right|}}, \ end{aligned}

(25)

where ({r}) is the total number of puncturing steps, i.e., (r = text{31}), and (i) is the position number of the output response, where (i=1,cdots ,10).

The procedure of needle insertion is recorded from the time of contact between the soft tissue and the needle tip to the time at which a depth of (90, mathrm{mm}) is reached inside the tissue. With the functional response Kriging model, the displacement in the (X) and (Y) directions can be predicted, and Figure 6 plots the displacements of N2 shown in Figure 2, which is close to the needle. In Figure 6, ‘Kriging value’ denotes the prediction results obtained using the functional response Kriging model, and ‘FEM’ denotes the results obtained with a finite element simulation. From Figure 6, the displacement fluctuations at the initial and final time steps are relatively large, while the fluctuations between these two time steps are smaller. The functional response Kriging model follows the law of tissue deformation during the needle insertion procedure. The model can also predict the displacements at other locations at different time steps very well. Figure 7 shows the displacements at N8, which is far away from the needle.

With Eq. (25), the mean residuals and relative mean residuals of the 10 locations are as listed in Table 2.

The fluctuation occurs at the initial time steps which is derived from the spring-back of the tissue surface as the needle punctured the tissue. And as the Kriging model is inherently a data-driven model which the interpolation function is applied to predict other unknown region, the prediction accuracy of the data which has large fluctuation would decrease. Overall, the Kriging model smooths the original finite element model, and the prediction accuracy of the initial stage is poor. The maximum residuals of the Kriging predictions and finite element calculations in the (X) and (Y) directions are (0.37, mathrm{mm}) and (1.24, mathrm{mm}), respectively, and both occur in the fifth time step, at which the needle tip punctures the surface of the soft tissue. The average residuals in the (X) and (Y) directions occur within (0.05, mathrm{mm}) and (0.12, mathrm{mm}), respectively, and the relative mean residual is at most (16%).

The runtime is also compared. Both the Kriging model and FEM model run on the same computer with macOS Mojave 10.14.6, a 2.7 Hz Intel Core i5 and 8 GB of RAM, and the program is written in MATLAB 2019b. The average runtime of 5 runs is recorded, with the average time for the Kriging-based method being 0.0294 s and that of the FEM model being 4.6912 s. Note that neither of the two reported times includes a visualization of the soft tissue deformation process. The reported results indicate that the Kriging model runs approximately 160 times faster than the FEM model.

### Adaptation of the Deformation Modeling to Step Variations

From the above analysis, the Kriging prediction model based on the observations at 31 time steps can reflect the soft tissue deformation at fixed (known) time steps. However, it is necessary to obtain the motion information of the sensitive area inside the soft tissue at any time step for both the path planning and the virtual training system. In this section, the adaptation to different time steps is studied.

Additional insertion time steps are expressed as ({d}_{1}=[1.0text{ mm},2.5text{ mm},5.0text{ mm},cdots ,90text{ mm}]), the insertion depth increment is 2.5 mm and the total depth achieved by the needle is (90, mathrm{mm}). The displacements of 10 observation locations at ({d}_{1}) insertion time steps are recorded to compare with the Kriging model. Forty sets of data randomly selected from 100 sets of records are used to construct the real-time Kriging prediction model of soft tissue deformation, and an independent set of data is used to test the model.

Figure 8 shows a comparison between the prediction and the computer-based experimental results at the N3 location. From Figure 8, the Kriging model with a functional response can fit the computer-based experimental results very well and can reflect the fluctuation trends at the initial and end time steps. The maximum residual occurs at the 6th insertion time step, and its maximum residuals and relative mean residuals in the (X) direction are (0.22text{ mm}) and (21%), respectively. The maximum residuals and relative residuals in the (Y) direction are (0.55text{ mm}) and (8%), respectively. The maximum residuals and relative mean residuals of the remaining nine observation positions are listed in Table 3.

From Table 3, the maximum residual in the (X) direction occurs at N1 and N7; its value is (0.31text{ mm}). The maximum residual in the (Y) direction appears at N2; its value is (0.92text{ mm}). The relative mean residual in the (Y) direction is less than (8%) and has a relatively high accuracy. For the real insertion procedure, the displacements in the (X) direction are very small, sometimes approaching zero; hence, the relative mean residual is larger than (20%). However, the absolute mean residual of each location is less than (0.31text{ mm}), which can meet the estimation accuracy.