Kinematic model and its parameter identification for cannula flexible needle insertion into soft tissue

Abstract

In minimally invasive surgery, flexible needle insertion is a popular application which has been extensively researched. However, needle steering is challenging for a bevel tip cannula flexible needle due to the nonholonomic constraints and the rebounds of the needle shaft when the needle tip is reoriented. We proposed a novel kinematic model for the bevel tip cannula flexible needle based on bicycle and unicycle models, taking consideration of the deflection of the bevel tip and the rebounds of the needle shaft. Aiming at different types of paths, forward kinematics of the model was analyzed and calculated. Each parameter of the kinematic models was identified based on the experimental data using the least square method. Furthermore, the changing rules of parameters were explored under different angles of the bevel tip. The experimental results show that the errors of the proposed kinematic models are within 2 mm, which is greatly reduced compared to the traditional bicycle or unicycle model, and fulfill a general clinical surgery.

Introduction

Minimally invasive surgery (MIS) is a brand new medical treatment, which has been developed rapidly in recent years. It wins good graces for its advantages of fewer traumas, rapider effectiveness and fewer complications, and so on. And it will be an inevitable developing trend of surgical treatments in modern clinical medicine. Needle insertion is a commonly used treatment in MIS, such as tissue biopsy and brachytherapy applied in prostate, liver, or brain. Nowadays, rigid steel needles are commonly used in the surgeries. However, due to the structural complexity of human body and uncertainty of practical operations, the traditional rigid needle can hardly meet the needs of complex surgeries, that is, it is not only difficult to rectify the deviated path (even though they often can be visualized), but also unable to steer clear of the obstacles (like bones) or sensitive tissues (like nerves and blood vessels). Targeting errors or failure arrivals to the target can diminish the effectiveness of therapies or bias the diagnosis results depending on the accuracy of needle placement. Thus, the needle steering technique is the key for the needle insertion procedures, while the kinematic model is the basic research for the steering, that is, we have to perform path planning and even control of the needle based on the kinematic model. ¹

Webster et al. ¹ proposed a bevel tip flexible needle (BTFN), which is usually made of nitinol (an alloy of nickel and titanium). The needle is supposed to be flexible enough to bend inside the tissue because of the lateral resistance acted on the bevel tip when inserted. Thereafter, other kinds of flexible needles with different types of bevel tips are generated.^2–4 Ko et al. ⁵ proposed a programmable bevel tip probe: the offset between two interlocked segments of the probe determines the bending direction and radius. Among all kinds of flexible needles, the BTFN proposed by Webster et al. has been widely accepted and extensively studied. However, the BTFN mainly has two drawbacks: one is that the curvature of the needle path is almost invariable in a certain tissue and the path can hardly be rectified if it is deviated; and the other is that the reorientation precision of the bevel tip is low inside the tissue, because there is a torsional friction along the needle shaft.

We have proposed a cannula flexible needle, which comprises a flexible cannula and a flexible stylet, as shown in Figure 1. The cannula and the stylet have relative translations and rotations. And both drawbacks above mentioned can be easily overcome by the cannula flexible needle: first, it generates a set of variable-curvature paths by inserting the stylet with different lengths of d out of the cannula; second, because the stylet and the tissue are separated by the cannula, the torsional friction along the stylet is greatly reduced; hence, the reorientation precision of the bevel tip is fundamentally improved.

Figure 1.

Cannula flexible needle.

Kinematic modeling for the BTFN steering is challenging due to the nonholonomic constraints of the needle and has been widely studied by scholars. DiMaio et al. formulated a needle Jacobian that described needle tip motion due to the needle base movement and proposed a model for the tissue based on a finite element method (FEM). However, the study was concerning a symmetric-tip needle. ⁶ Alterovitz et al. simulated bending of the BTFN based on FEM. The simulation modeled the effects of the needle tip and frictional forces on soft tissues defined by a two-dimensional mesh. ⁷ Chentanez et al. ⁸ proposed an interactive needle insertion simulator based on FEM that modeled the coupling between the BTFN and deformable tissue. However, the FEM is time-consuming and cannot generate any general analytic equations. Misra et al. ⁹ first developed a mechanics-based model for the BTFN steering, which can be used to predict the needle behavior based on mechanical and geometrical properties of the needle. Wang et al. ¹⁰ also developed a mechanics-based model to predict the needle deflection, utilizing the beam theory and the Rayleigh–Ritz method. Datla et al. ¹¹ presented an energy-based model to predict needle deflection of an active bevel-tipped needle when inserted into tissue. Khadem et al. ¹² proposed mechanics-based model for simulating a planar needle insertion in soft tissue based on the beam theory. Jiang and Wang ¹³ also proposed a mechanics-based model for the BTFN insertion based on the quasi-static principle and the dynamics of the needle motion. However, the mechanics-based or the energy-based models are complex and also time-consuming, which may not be suitable for a real-time prediction or path planning.

Webster et al. ¹ proposed a bicycle model (BM) and a unicycle model (UM) for the BTFN, considering the path of the needle was equivalent to a path of a bicycle or a unicycle. Other kinematic studies on the BTFN were similar to the BM or the UM.^14–16 These kinematic models are the analytic equations, which are appropriate for calculation of paths. Therefore, a great deal of path planning has been studied based on the BM or the UM.^17–19 However, these models neglect the elasticity of the needle shaft, assuming that the needle is completely flexible and the needle tip “stays in place” during reorientation. As a matter of fact, the needle tip will rebound to some extent during reorientation, which definitely causes a large error between the model and the actual path. Moreover, different types of conjunctions between segments of a path, for example, arc-arc (A-A), line-arc (L-A), or arc-line (A-L), will have different rebounds. Although literatures^1,10,20 proposed models to predict double bendings (i.e. A-A), they did not take account of the rebound in their models. In our previous work, we proposed a deflection unicycle model with rebound (UMR), which was a modified UM, and calculated the kinematics of a path with A-A type for the BTFN. ²¹ However, to our best knowledge, there is no existing study on kinematic models sorted by different types of paths, and no previous work on model parameter identifications considering the rebounds of L-A or A-L for the BTFN. Especially, there is no previous work discussing the changing rules of each parameter. In this article, we will propose a novel kinematic model for the cannula flexible needle. We are going to establish the models for the different types of paths (including the single arc, A-A, L-A, and A-L) considering the deflection of the needle tip and the rebounds of the needle shaft, calculate the kinematics, identify the model parameters, and even explore the changing rules of parameters for the kinematic models of the cannula flexible needle.

Kinematic model and calculation of cannula flexible needle

Establishment of kinematic model

Motion analysis

Different with the traditional BTFN, which has 2 degrees of freedom (DOF)s, that is, insertion and rotation, there are 3 DOFs for the cannula flexible needle: translation u₁ for the cannula, translation u₂ and rotation u₃ for the stylet. The radius of the needle path is varied by the length d inserted out of the cannula. The larger the length d gets, the smaller the radius will be. Therefore, the relative movement of u₁ and u₂ changes the curvatures of paths, while u₃ determines the directions of bending. However, after adjusting the length d, we usually input u₁ and u₂ simultaneously in the same velocity in order to get a certain radius of path; thus, we denote u₁ and u₂ as u₁₂. If u₁₂ and u₃ are input separately, an arc-based path will be generated, see Figure 2(a). If u₁₂ and u₃ are input simultaneously, and the speed of u₃ is much smaller than u₁₂, a helix-based path will be generated, see Figure 2(b); however, if the speed of u₃ is much larger than u₁₂, a linear path will be generated, see Figure 2(c). ²²

Figure 2.

Motion forms: (a) arc-based path, (b) helix-based path, and (c) linear path.

Modified unicycle model

First, let us recall the BM and the UM as a reference (see Figure 3). The functions of the needle path for the BM and the UM are as follows

$$\{\begin{array}{c}\mathit{p}=f(\phi ,l_1,l_2)\mathrm{for}\mathrm{BM}\hfill \\ \multicolumn{1}{c}{\mathit{p}=f(r)\mathrm{for}\mathrm{UM}}\end{array}$$

(1)

Figure 3.

Bicycle and unicycle models: (a) bicycle model and (b) unicycle model.

The experiments showed that the BM was more accurate than the UM for a single arc. ¹ In order to simplify the kinematic calculation and not to diminish the accuracy of the kinematic model, we propose a new modified unicycle model (MUM), which is similar to BM for a single-arc path. However, in contrast with the BM, the number of parameters in the MUM is reduced into two

$$\mathit{p}=f(r,\phi )$$

(2)

Because of the gap between the needle and its support, when the cannula flexible needle is bending, the bevel (i.e. the wheel) will be deflected at the end of the support. Thus, the arc path of the needle is not the original unicycle path; instead, it is virtually an arc with a deflection of the wheel (rotation about x-axis by φ), as shown in Figure 4. Therefore, we should modify it with a deflection φ. The MUM is only suitable for the path start with an arc (A-A and A-L), but not suitable for that start with a line (L-A).

Figure 4.

Modified unicycle model. The bevel is deflected by an angle φ with respect to the needle axis.

Models with rebounds

Although the MUM (or the BM) is more accurate than the UM, it remains intolerant errors for non-single-arc paths. The main cause is under the assumption of “staying in place” of the needle tip during reorientation. As a matter of fact, when the needle is rotated inside the tissue, the bevel direction changes, and the effective interaction area between the bevel and the tissue will be reduced; thus, the resistance acted on the tip is also reduced. This will lead to a certain rebound of the needle tip, as shown in Figure 5. Hence, a large error occurs between the actual path and the prediction of the MUM, as shown in Figure 6. Thus, we have to take account of the rebound after reorientation. Moreover, different types of paths will have the different rebounds. So we propose a UMR and a modified unicycle model with rebound (MUMR). The superiority of the modification of the rebound is obvious, and the experimental contrasts are available in Zhang and Zhao. ²¹

Figure 5.

Rebound of the needle: (a) before rotation and (b) after rotation.

Figure 6.

Comparison of actual path and prediction of MUM.

Of course, different rotation angles of the needle shaft will result in the different rebounds, which is more complicated; however, in this article, we focus on the study of the path in a plane (rotation with 180°).

Now, we sort the paths into three types, that is, A-A, L-A, and A-L, as shown in Figure 7. The needle tip is inserted to the point N. After rotation, the segment CN of the original path will rebound to be the segment CD, as shown in Figure 7(a) and (b). While for L-A, the needle is inserted to the point N, after insertion of the arc, the segment CN will be pulled into CD because of the lack supporting of the tissue. So we define CDs as rebounds (noted as b ). The rebounds should connect the conjunctive segments of a path smoothly, that is, the segment CD is tangent to both arcs in A-A (noted as b₁, see Figure 7(a)), or is tangent to the first segment of arc and identical with the following line in A-L (noted as b₂, see Figure 7(b)), or is tangent to the first segment of line and identical with the following arc in L-A (noted as b₃, see Figure 7(c)), respectively. Thus, we have to take account of the rebounds b  = (b₁, b₂, b₃) into the UM and the MUM. Then the UMR and the MUMR are obtained. And their path functions are as follows

$$\mathit{p}=\{\begin{array}{c}f(r,\mathit{b})\mathrm{for}\mathrm{UMR}\hfill \\ \multicolumn{1}{c}{f(r,\phi ,\mathit{b})\mathrm{for}\mathrm{DUMR}}\end{array}$$

(3)

Figure 7.

Rebounds of different types of paths: (a) A-A typed path with rebound, (b) A-L typed path with rebound, and (c) L-A typed path with rebound.

Kinematic models for different types of paths

To sum up, the kinematic models of the cannula flexible needle are the modified UMs, which are improved with the deflection and the rebounds, that is, if the first segment of the path is an arc (including a single arc), the path needs to be modified with the deflection to the UM (i.e. MUM); however, it need not for the path that the first segment is a line (i.e. UM). If the path consists of multiple segments, it needs to be modified with the rebounds (i.e. UMR or MUMR). And the general function of the kinematic model is formulated as follows

$$\mathit{p}=f(r,\phi ,\mathit{b})$$

(4)

where r is the radius of a segment; φ is the deflection angle of the wheel, if the first segment of path is a line, φ = 0, and the model is degraded into the UMR; b is the rebound, if the path consists of only one segment, b  = 0, and the model is degraded into the MUM. Detail relationships between paths and models are as follows

$$\begin{array}{c}\mathrm{Kinematic}\hfill \\ \multicolumn{1}{c}{\mathrm{Models}}\end{array}\{\begin{array}{c}\begin{array}{c}\mathrm{Single}\hfill \\ \multicolumn{1}{c}{\mathrm{segment}}\end{array}\{\begin{array}{c}\mathrm{Line}:\mathrm{UM}\hfill \\ \multicolumn{1}{c}{\mathrm{Arc}:\mathrm{MUM}}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}\mathrm{Multiple}\hfill \\ \multicolumn{1}{c}{\mathrm{segments}}\end{array}\{\begin{array}{c}\mathrm{The}\mathrm{first}\mathrm{segment}\mathrm{is}\mathrm{a}\mathrm{line}:\mathrm{UMR}\hfill \\ \multicolumn{1}{c}{\mathrm{The}\mathrm{first}\mathrm{segment}\mathrm{is}\mathrm{an}\mathrm{arc}:\mathrm{MUMR}}\end{array}}\end{array}$$

Kinematic calculation

Calculation for unicycle kinematic model

As shown in Figure 3(b), the configuration of the needle tip frame B relative to the world frame A can be described by a homogeneous transformation matrix

$${\mathit{g}}_{\mathrm{ab}}=\left[\begin{array}{cc}{\mathit{R}}_{\mathrm{ab}}\hfill & {\mathit{p}}_{\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]\in \mathrm{SE}(3)$$

(5)

where R _ab  ∈ SO(3) is the rotation matrix, and ${\mathit{p}}_{\mathrm{ab}}\in {\mathbb{R}}^3$ is the position vector of frame B relative to frame A.

In frame B, the twist $\hat{\xi }\in \mathrm{se}(3)$ is formulated as

$$\hat{\xi }={\left[\begin{array}{c}\mathit{v}\hfill \\ \multicolumn{1}{c}{\mathbf{\omega }}\end{array}\right]}^{\wedge }=\left[\begin{array}{cccc}0\hfill & -u_3\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{u_3}& 0\hfill & -u_{12}/r\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& u_{12}/r\hfill & 0\hfill & u_{12}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(6)

where v is the instantaneous linear velocity of frame B relative to frame A, $\mathit{v}=[\begin{array}{ccc}0\hfill & 0\hfill & {u_1}_2\hfill \end{array}{]}^{\mathrm{T}}$ ; ω is the angular velocity, $\mathbf{\omega }\mathbf{=}\mathbf{[}\begin{array}{ccc}{\mathbf{u}}_{12}\mathbf{/}\mathbf{r}\hfill & 0\hfill & {\mathbf{u}}_3\hfill \end{array}{\mathbf{]}}^{\mathrm{T}}$ ; r is the radius of a segment; and ^ is the wedge operator to form a matrix into se(3) out of a given vector in ${\mathbb{R}}^6$ .

Then, the homogeneous transformation matrix can be rewritten into an exponential form

$${\mathit{g}}_{\mathrm{ab}}(t)={\mathit{g}}_{\mathrm{ab}}(0)\exp (\hat{\xi }t)$$

(7)

where g _ab (0) is the initial configuration of the needle before insertion (frame B relative to frame A); t is the insertion time. Details could be attained in Murray et al. ²³

Kinematic calculation for single segment

For a single segment, all we have to consider is to modify with the deflection and do not need to consider the rebounds. According to equation (7), the homogeneous transformation for the MUM can be obtained as follows

$${\mathit{g}}_{\mathrm{ab}}(t)={\mathit{g}}_{\mathrm{ab}}(0){\mathit{g}}_0\exp (\hat{\xi }t)$$

(8)

where g _ab (0) is the initial configuration of the needle; g ₀ is the modified needle configuration with regard to the UM, if the single (or the first) segment is an arc, when t > 0, g ₀ is the transformation matrix of rotation about x-axis by a deflection angle φ

$${\mathit{g}}_0={\mathit{g}}_x(\phi )=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cos \phi \hfill & -\sin \phi \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \sin \phi \hfill & \cos \phi \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(9)

However, when t = 0, ${\mathit{g}}_0={\mathit{I}}_4$ , where I ₄ is a unit matrix with 4 by 4, that is, the needle is in the initial configuration.

While if the single (or the first) segment is a line, it need not be modified with the deflection; thus, ${\mathit{g}}_0={\mathit{I}}_4$ . Actually, it is the UM now. And the term $\exp (\hat{\xi }t)$ in equation (8) will be approximate to

$$\exp (\hat{\xi }t)\approx {\mathit{g}}_l(L){\mathit{g}}_z(\gamma )$$

(10)

where ${\mathit{g}}_l(L)$ is the transformation matrix of translation with L along z-axis, L = vt; ${\mathit{g}}_z(\gamma )$ is the transformation matrix of rotation by $\gamma$ about z-axis, $\gamma$ is the rotation angle of the needle shaft, $\gamma =\omega t-2k\pi \in [\begin{array}{cc}0\hfill & 2\pi )\hfill \end{array}$ , k = 0, 1, 2, ….

Kinematic calculation for multiple segments

Kinematic model for multiple segments is rather complicated. Not only the deflection but also the rebounds have to be taken into account. Since the whole path can be discretized into several segments, let ${\mathit{g}}_i=\exp ({\hat{\xi }}_i{t}_i)$ be the transformation matrix for the ith segment, where ${\hat{\xi }}_i$ and t_i are the twist and the insertion time for the ith segment, respectively.

If the rebounds are neglected, the forward kinematics after insertion of n segments is formulated as

$${\mathit{g}}_{\mathrm{ab}}(T)={\mathit{g}}_{\mathrm{ab}}(0){\mathit{g}}_0{\mathit{g}}_1{\mathit{g}}_2\dots {\mathit{g}}_n$$

(11)

where T is the total insertion time of the whole path, $T=\sum _{i=1}^n{t}_i$ ; g _ab (0) is the initial configuration; g ₀ is the modified configuration of the needle.

As a matter of fact, there are rebounds existing in the path; therefore, we have to modify each segment until the last one, that is, to add transformations of the rebounds between segments. As shown in Figure 7, let BN be the ith segment with the transformation matrix g _i . There is a rebound b_j (j is to sort the type of rebounds, j = 1, 2, 3) between the ith and the i + 1th segments. Thus, the insertion movement of the two segments is equivalent to this procedure: first to insert the needle from the point B to the point C for the ith segment, then to experience the rebound (moving from C to D); then to reorientate and to insert the i + 1th segment. Therefore, the insertion of the ith segment has to be modified as

$${\tilde{\mathit{g}}}_i=\exp ({\hat{\xi }}_i(t_i-t_b))$$

(12)

where t_b is the modified time, that is, the insertion time for the length of the rebound, and different rebounds possess the different modified times, $t_b=b_j/u_{12}$ , j = 1, 2, 3.

The transformation for the rebounds is as follows

$${\mathit{g}}_b=\{\begin{array}{c}{\mathit{g}}_l(b_j){\mathit{g}}_z(\gamma \prime )\mathrm{for}\mathrm{A}\mathrm{-}\mathrm{A}\mathrm{or}\mathrm{A}\mathrm{-}\mathrm{L}\hfill \\ \multicolumn{1}{c}{{\mathit{g}}_z(\gamma \prime ){\mathit{g}}_r(t_b)\mathrm{for}\mathrm{L}\mathrm{-}\mathrm{A}}\end{array}$$

(13)

where ${\mathit{g}}_z(\gamma \prime )$ is the transformation matrix of rotation about z-axis, $\gamma \prime$ is the rotation angle of the needle shaft before insertion for the next segment (it is 180° for a planar path); ${\mathit{g}}_r(t_b)$ is the transformation matrix of an arc, which can be obtained by equation (7) when the input u₃ = 0 in the twist, t_b is the modified time.

Thus, the kinematics of the MUMR is as follows

$${\mathit{g}}_{\mathrm{ab}}={\mathit{g}}_{\mathrm{ab}}(0){\mathit{g}}_0{\tilde{\mathit{g}}}_1{\mathit{g}}_b{\tilde{\mathit{g}}}_2{\mathit{g}}_b\dots {\tilde{\mathit{g}}}_k{\mathit{g}}_b\dots {\mathit{g}}_n$$

(14)

Parameter identification for kinematic models

Parameter identification in single arc

From equation (4), we can conclude that the center of the arc is a function with the parameters of r and φ

$$\mathit{c}=f(r,\phi )$$

(15)

where r is the radius of the arc; φ is the deflection angle.

In order to identify the two parameters, we use the MUM to fit the single arc. Let the starting point be p₀ (0, y₀, z₀), and a sample point in the path be p_i (0, y_i, z_i), as shown in Figure 8. We assume that the needle axis is initially aligned with the z-axis of frame A and the center of the arc O (x_c, y_c, z_c) is formulated as follows

$$\mathit{c}={\left[\begin{array}{ccc}{x}_c\hfill & y_c\hfill & z_c\hfill \end{array}\right]}^T={\left[\begin{array}{ccc}0\hfill & y_0-r\cos \phi \hfill & z_0-r\sin \phi \hfill \end{array}\right]}^T$$

(16)

Figure 8.

Position relations in single arc.

Since the position of p_i is in the arc centered at O; thus, we have

$$z_i=f_c(y_i)=\pm \sqrt{r^2-{(y_i-y_c)}^2}+z_c$$

(17)

where “±” depends on which side of p_i located at O; if z_i > z_c, it is “+”; on the contrary, it is “–.”

The least square method is adopted, and the both parameters of r and φ can be identified by

$$F=\min \sum _{i=1}^n{({\tilde{z}}_i-f_c({\tilde{y}}_i))}^2$$

(18)

where ${\tilde{y}}_i$ and ${\tilde{z}}_i$ are the experimental coordinates of p_i.

Although the calculation in equation (17) includes “±” signs, the sign is certain with the relative position of p_i to O; therefore, the calculation is definite. Thus, the identification of the parameters is unique, and the same as the following parameter identifications.

Parameter identification in A-A

We use the MUMR to fit the path A-A in order to identify the parameters of r, φ, and b₁. The center O₁ of the first arc is identical with that of the single arc, that is

$$c_1=[\begin{array}{ccc}0\hfill & y_{c_1}\hfill & z_{c_1}\hfill \end{array}{]}^{\mathrm{T}}=c$$

(19)

Let the length of the first segment BN be S₁, and the second segment be S₂, which are assumed to be known in advance according to the insertion length (see Figure 7(a)), so we have

$${\alpha }_1=\frac{S_1-b_1}{r}$$

(20)

where S₁ is the length of the first segment; b₁ is the rebound of A-A; r is the radius of the arcs.

In the triangle O₁CE, we have

$${\alpha }_2=\arctan \left(\frac{b_1}{2r}\right)$$

(21)

So

$$\beta =\phi +{\alpha }_1+{\alpha }_2$$

(22)

where φ is the deflection angle.

The distance between the two arc centers

$$l=2\sqrt{{\left(\frac{b_1}{2}\right)}^2+r^2}$$

(23)

Thus, the coordinates of O₂ can be expressed by the position of O₁

$${\mathit{c}}_2={\left[\begin{array}{ccc}0\hfill & y_{c_2}\hfill & z_{c_2}\hfill \end{array}\right]}^T={\left[\begin{array}{ccc}0\hfill & y_{c_1}+l\cos \beta \hfill & z_{c_1}+l\sin \beta \hfill \end{array}\right]}^T$$

(24)

For a point p_i(0, y_i, z_i) in the path, the relationship between y_i and z_i is as follows

$$z_i=f_{c_j}(y_i)=\pm \sqrt{r^2-{(y_i-y_{c_j})}^2}+z_{c_j}$$

(25)

where j represents the jth segment, j = 1, 2; “±” depends on which side of p_i located at O_j, for the first segment, $p_i\in S_1$ , if $z_i>z_{c_1}$ , it is “+”; on the contrary, it is “–”; and similarly for the second segment, $p_i\in S_2$ , if $z_i>z_{c_2}$ , it is “+”; on the contrary, it is “–.”

Then, the least square method is adopted again. By substituting equation (25) into equation (18), the three parameters of r, φ, and b₁ can be identified.

In the identification, although the path consists of two segments, the coordinate data are definite for the equations of each segment. Thus, the identification of the parameters is unique also, and the same as the following parameter identifications.

Parameter identification in A-L

We also use the MUMR to fit the path A-L in order to identify the parameters of r, φ, and b₂. The center O of the first arc is also identical with that of the single arc, as formulated in equation (19).

For a point p_i(0, y_i, z_i) in the first segment, the relationship between y_i and z_i is the same as equation (17). Also let the length of the first segment BN be S₁ (see Figure 7(b)), so we have

$$\alpha =\frac{S_1-b_2}{r}$$

(26)

where S₁ is the length of the first segment; b₂ is the rebound of A-L; r is the radius of the arc

$$\beta =\phi +\alpha$$

(27)

where φ is the deflection angle.

Thus, the position of the point C is formulated as follows

$${\mathit{p}}_c={\left[\begin{array}{ccc}{p}_x\hfill & p_y\hfill & p_z\hfill \end{array}\right]}^T={\left[\begin{array}{ccc}0\hfill & y_c+r\cos \beta \hfill & z_c+r\sin \beta \hfill \end{array}\right]}^T$$

(28)

The slope of the line is

$$k=-\tan (\beta )$$

(29)

And, for a point p_i(0, y_i, z_i) in the linear segment, the relationship between y_i and z_i is as follows

$$y_i=k(z_i-p_z)+p_y$$

(30)

Then, by substituting equations (17) and (30) into equation (18) for each segment, the three parameters of r, φ, and b₂ can be identified.

Parameter identification in L-A

We use the UMR to fit the path L-A in order to identify the parameters of r, φ, and b₃, as shown in Figure 7(c). Assuming that the needle is inserted in horizon, the first linear segment is formulated as

$$y=y_0$$

(31)

And the center O of the arc (the second segment) is formulated as follows

$${\mathit{c}}_2={\left[\begin{array}{ccc}0\hfill & y_{c_2}\hfill & z_{c_2}\hfill \end{array}\right]}^T={\left[\begin{array}{ccc}0\hfill & y_0-r\hfill & S_1-b_3\hfill \end{array}\right]}^T$$

(32)

where S₁ is the length of the first segment; b₃ is the rebound of L-A; r is the radius of the arc.

And the relationship between y_i and z_i is the same as equation (17). Then, by substituting equations (17) and (31) into equation (18) for each segment, the three parameters of r, φ, and b₃ can be identified.

Experimentation and result

Parameter identification

The needle experimental system consists of an insertion mechanism, a control system, and a block of phantom tissue, as shown in Figure 9. The insertion mechanism is composed of two screw–nut mechanisms, which are overlapped. The bottom one is to realize the insertion of both the cannula and the stylet by motor I. The upper one is to realize the insertion of the stylet relative to the cannula (adjusting the length d of the stylet out of the cannula) by motor II. And motor III which is mounted on the upper screw-nut drives the needle shaft spin.

Figure 9.

Experiment setup of the cannula flexible needle.

The cannula flexible needle used in the experiment was as follows: the stylet was a cylinder nitinol wire with the diameter of 1 mm, and with the bevel tip of 30°; and the cannula was a tube with the inner diameter of 1.1 mm and the outer diameter of 1.3 mm. In the experiment, the length d was set to zero, that is, only the bevel tip remains out of the cannula. And the phantom tissue was adopted as agar.

Although the parameter values (including r, φ, and b ) of the kinematic models will be affected by the diameter of the needle, the length of d, the angle of the bevel tip, and even the properties of the needle or the tissue, they are constant for a specific situation, the paths are similar, and the considerations of kinematic modeling are the same.

The velocities in the experiment were defined as v = 0.009 m/s and ω = 3.5 rad/s, respectively. The insertion lengths for each type of paths were as follows: 110 mm for the single arc, 50 mm for the first segment, and 60 mm for the second, respectively. Rotation angle for reorientation of A-A was 180° for sake of the planar paths. Each type of path was inserted for five times. One of the insertions for each type is shown in Figure 10, which was acquired by a digital camera. And the trajectories were obtained by using an image processing.

Figure 10.

Experiments for parameter identifications: (a) single arc, (b) A-A, (c) A-L, and (d) L-A.

Assuming that all the radii of arcs in all types of paths remain the same, and in order to obtain unanimous parameters, we put all the experimental data together and fit all the parameters at the same time for the parameter identifications. The results are shown in Figure 11 and Table 1.

Figure 11.

Fitting results: (a) single arc, (b) A-A, (c) A-L, and (d) L-A.

Table 1.

Parameter identification results.

r (mm)	φ (°)	b 1 (mm)	b 2 (mm)	b 3 (mm)	Maximum error	RMS of errors
63.7	2.1	12.1	5.9	4.0	1.7223	0.5471

RMS: root mean square.

From the results, we can conclude that the kinematic models fit the actual paths very well, with a little maximum error (less than 2 mm) and a little root mean square (RMS) of errors. This proves the validity of the proposed kinematic models with the presented parameters.

Kinematic model comparison

In order to compare the proposed models with the traditional bicycle or unicycle models, and prove the validity of the proposed parameters (the deflection and the rebounds), we also did experiments for the comparison. MUMR is with both the deflection and the rebounds, UMR is with the rebound but without the deflection, MUM is with the deflection but without the rebounds, and UM is without any of the parameters.

We chose the single-arc path and the A-A typed path, because they possess both parameters in the proposed models for the comparison. The results are shown in Figure 12 and Table 2.

Figure 12.

Comparison of the models: (a) MUMR, (b) UMR, (c) MUM, and (d) UM.

Table 2.

Results of the models.

Models	RMS of errors	Maximum error
MUMR	0.5480	1.5882
UMR	0.9527	2.1331
MUM	2.1849	6.3718
UM	2.4616	8.3306

RMS: root mean square; MUMR: modified unicycle model with rebound; UMR: unicycle model with rebound; MUM: modified unicycle model; UM: unicycle model.

From the results, we can conclude that the models with rebound are much superior to those without it, and the models with deflection are superior to those without it. This proves the rationality of the proposed parameters.

Changing rules of parameters

As above noted, the parameters’ values of the models will vary under different bevel angles of the needle tip. In order to discuss the changing rule of each parameter under different bevel angles, we perform experiments with eight different angles of the needle tip, that is, 10°, 20°, 30°, 40°, 50°, 60°, 70°, and 80°.

Rule of radius r

The change of the radius r is continuous, and it seems to be symmetrical and with a minimum value, which accords with the change of cosecant function very much, as shown in Figure 13. Thus, the regression result of r is as follows

$$r=28.2716\csc 2\theta +32.4635,0°<\theta ⩽90°$$

(33)

Figure 13.

Cosecant regression for r.

The maximum error and the RMS of errors are 1.6538 and 1.0610 mm, respectively. This proves the rationality of the regression rule for r.

Rule of deflection φ

We respectively adopt a quadratic curve and a sine curve to regress the experimental data under different bevel angles, and the results are as follows

$$\{\begin{array}{c}\phi =-0.0014{\theta }^2+0.1254\theta -0.4946,\mathrm{quadratic}\mathrm{curve}\hfill \\ \multicolumn{1}{c}{\phi =2.3097\sin 2\theta ,\mathrm{sine}\mathrm{curve}}\end{array}$$

(34)

As shown in Figure 14 and Table 3, both regressions fit the data well, but they have the superiority of their own, the error of the quadratic result is a little less than that of the sine result; however, the sine result accords with the actual limitations for 0° and 90°, that is, theoretically, there will be no deflection (the value should be zero) if the bevel angle is 0° or 90°.

Figure 14.

Regression results of φ under different θ.

Table 3.

Regression errors of φ.

Regression	Maximum error (mm)	RMS of errors (mm)
Quadratic curve	0.3589	0.1717
Sine curve	0.3846	0.1816

RMS: root mean square.

Rule of rebounds b

We adopt quadratic curves to regress the experimental data under different bevel angels, and the regression results of b (b₁, b₂, and b₃) are as follows

$$\{\begin{array}{c}{b}_1=0.0070{\theta }^2-0.6348\theta +25.0554\hfill \\ \multicolumn{1}{c}{b_2=0.0037{\theta }^2-0.3299\theta +11.6839}\\ \multicolumn{1}{c}{b_3=0.0026{\theta }^2-0.2262\theta +8.2464}\end{array}$$

(35)

As shown in Figure 15 and Table 4, the results show that the regressions fit the data well. This proves the rationality of the regression rules for b .

Figure 15.

Regression results of b under different θ.

Table 4.

Regression errors of b .

Rebounds	Maximum error	RMS of errors
b 1	0.4530	0.2601
b 2	0.3994	0.2325
b 3	0.2679	0.1625

RMS: root mean square.

Conclusion and future work

Based on the BM and the UM, we proposed the modified kinematic models for the cannula flexible needle, taking consideration of the deflection of the bevel and the rebounds of the needle shaft when reorientated. Then, the kinematic models are sorted according to the different types of paths. Forward kinematics of each model is analyzed and calculated. The identifications for each parameter (arc radius r, deflection angle φ, and rebounds b ) of the models are performed based on the experimental data using the least square method. Furthermore, the changing rules of the parameters under different bevel angles are discussed using a regression method. The experimental results prove that the proposed models are appropriate for the cannula flexible needle, and the parameter identifications are rational.

Moreover, the proposed models are suitable for any planar path no matter how many segments it possesses. Because the path types mentioned in the manuscript are the basic elements, and a path is the conjunction of the basic ones. In addition, the model parameters remain the same, and the kinematics can be calculated, accordingly.

Although the real biological tissues are inhomogeneous and exhibit non-linear properties which are quite different from the agar phantoms adopted in our experiments, the kinematic models are similar, merely with the different values of parameters. Although this may cause the larger errors of the kinematic models, we could estimate and modify the parameters in a real-time control with a medical imaging technology and use a dynamic path planning algorithm to rectify the deviated path and diminish the errors.

In future work, we will study the sensitivity of the model with respect to these parameters, explore the relationships between the model parameters and the properties of needle or tissue, and further explore the rebounds with any other angles (not limited to 180°). Moreover, the study of path planning based on the proposed kinematic models for the needle steering is also valuable work.