Abstract
A variety of medical treatment and diagnostic procedures rely on flexible instruments such as catheters and endoscopes to navigate through tortuous and soft anatomies like the vasculature. Knowledge of the interaction forces between these flexible instruments and patient anatomy is extremely valuable. This can aid interventionalists in having improved awareness and decision-making abilities, efficient navigation, and increased procedural safety. In many applications, force interactions are inherently distributed. While knowledge of their locations and magnitudes is highly important, retrieving this information from instruments with conventional dimensions is far from trivial. Robust and reliable methods have not yet been found for this purpose. In this work, we present two new approaches to estimate the location, magnitude, and number of external point and distributed forces applied to flexible and elastic instrument bodies. Both methods employ the knowledge of the instrument’s curvature profile. The former is based on piecewise polynomial-based curvature segmentation, whereas the latter on model-based parameter estimation. The proposed methods make use of Cosserat rod theory to model the instrument and provide force estimates at rates over 30 Hz. Experiments on a Nitinol rod embedded with a multi-core fiber, inscribed with fiber Bragg gratings, illustrate the feasibility of the proposed methods with mean force error reaching 7.3% of the maximum applied force, for the point load case. Furthermore, simulations of a rod subjected to two distributed loads with varying magnitudes and locations show a mean force estimation error of 1.6% of the maximum applied force.
1 Introduction
Nowadays, flexible instruments and robots are increasingly being used for medical treatment and diagnostic purposes. Thanks to their compliant nature they can navigate through tortuous paths and adjust their shape according to the surrounding environment. Examples of such medical procedures are endovascular aneurysm repair, angioplasty and stenting, thrombolytic therapy, endoscopic gastroscopy, and colonoscopy ; ; . All these procedures require navigating a flexible instrument through a deformable lumen. For example, in endovascular aneurysm repair a stent catheter is moved through the vaculature. Exerting high forces on soft tissue, for short or extended periods of time, may cause a diversity of complications including perforation, symptomatic intracerebral hemorrhage, bleeding and ischemic complications ; ; . As a consequence, the demand for sensorized instruments relaying crucial real-time information such as: position, shape, and interaction force is ever increasing ; ; ; ; .
The problem of estimating externally applied forces on robotic manipulators is not new. There has been an abundance of works carried out on rigid link robots (; ; ; ), and more recently, on soft/continuum robots (; ; ; ). Regarding the latter, most of the previous works focused on estimating the contact force at the tip of the manipulator/instrument, rather than along the whole body. For example, Rucker et al. employed an Extended Kalman Filter (EKF) approach and the Cosserat rod model to estimate tip forces using pose measurements and a kinematic-static model of the robot (). Similarly, Hasanzadeh et al. made use of pose measurements, and introduced their own quasi-static piecewise circular arc model of the catheter to estimate the tip force (Hasanzadeh and Janabi-Sharifi, 2016). Back et al. use shape information and a simplified Cosserat rod model which can be rapidly solved using an iterative optimization algorithm to estimate forces at the tip of a catheter (). Hooshiar et al. estimate tip forces by employing Bezier spline shape approximations and an inverse Cosserat rod model (). Sadati et al. employed force sensor readings at the base of a continuum appendage and the Cosserat rod model for tip load estimation (). Further examples of previous works on tip force estimation using kinematic models can be found in (; ).
As tip force estimation was well understood in recent years, increasing attention was directed towards estimating forces along the length of the flexible instrument body. Qiao et al. made use of curvature measurements using multi-core fibers (MCFs) inscribed with fiber Bragg grating (FBG) sensors, a Cosserat rod model, and linear segmentation of curvature to estimate 2D lateral point forces subjected to an elastic rod body. The feasibility was tested on a Nitinol rod with two point loads (). In a later work, Qiao et al. employed an extended Kalman filter and pose information in combination with a Cosserat rod model to estimate point loads (). The approach was validated using two point loads on a thin Nitinol rod. Aloi et al. proposed a method to estimate distributed loads on an elastic rod using a Cosserat rod model and constrained nonlinear optimization. The method was tested using one, two and three point loads in diverse lateral directions along the cross-sectional plane. The paper focused on the feasibility of the proposed method but had a large computational cost and limited use in real-time applications (). Heunis et al. () and Venkiteswaran et al. () modelled a continuum manipulator subjected to multiple external loads using a pseudo rigid body (PRB) model. The location of external loads is assumed to be known a priori. Heunis et al. made use of three-dimensional shape information provided by FBG sensors inside the catheter. Ryu et al. employed shape information, a redundant kinematics model, and an extension of the Cosserat rod model to estimate point loads (). Lastly, Massari et al. combined machine learning and FBG sensors to develop a tactile sensor capable of estimating the magnitude and location of normal loads exerted on the surface (). The sensor, in theory, can be attached onto the external surface of flexible instruments. Although promising, it is limited with regards to the area of which a force can be exerted upon. In other words, it cannot cover the full circumferential surface of a circular instrument.
In this paper, two new approaches to estimate the number, location, and magnitude of external forces applied onto a flexible and elastic instrument body are proposed. Both approaches rely on Cosserat rod theory as the underlying modelling principle. Additionally, both approaches depart from the availability of curvature measurements. In this work, an MCF inscribed with FBGs is embedded within the instrument’s central working channel to obtain discretized curvature measurements along its length. However, any other method or sensing technology that could also provide similar curvature measurements would suffice. The first estimation approach makes use of piecewise polynomial curvature segmentation, while the second approach employs model-based optimization. The strategies are: 1) initially developed and experimentally tested for two-dimensional point forces applied anywhere along the instrument’s body, and 2) further adapted to estimate both point forces and distributed forces, where a simulation study is provided to prove its feasibility. Note that neither prior knowledge or estimate is required regarding the location of the forces nor an initial identification of whether it is a point force or a distributed force.
According to the above, the contributions of this paper provide the following:
1. force estimation methods that require no prior knowledge of the force interaction, but only knowledge of the two-dimensional curvature profile
2. three new piecewise polynomial-based curvature segmentation algorithms to estimate the number, location, and magnitude of external point forces to high accuracy
3. two new estimation-based approaches that solely employ curvature measurements to estimate the number, location, and magnitude of external point and distributed forces to high accuracy, which can be readily extended to include any types of loads
4. static and quasi-static experimental and simulation-based validations of the proposed methods up to three loads applied simultaneously
The rest of this paper is organized as follows: Section 2 briefly introduces the principles of curvature measurement using FBGs and the basic equations of Cosserat rod theory that are used in this work. Section 3 describes the two force estimation approaches in detail. Section 4 covers the experimental setup and the various experimental test configurations. Section 5 presents the results that were obtained. Section 6 will cover the method used for distributed force estimation and a corresponding simulation study. Finally, section 7 provides a conclusion and discussion of the proposed estimation approaches, and comments on future work prospects.
2 Sensing and Model Principles
2.1 Curvature Sensing
FBGs detect variation of strain based on the change of periodicity of a grating, and its refractive index. The Bragg wavelength is the wavelength of the light that is reflected back from the grating. The total change in strain can be due to a mechanical strain , or thermal expansion due to a temperature change . Considering small temperature shifts, the change in Bragg wavelength can be expressed as:where is the grating’s unstrained Bragg wavelength, is the change in strain with respect to the unstrained state, and and are the strain and temperature sensitivity coefficients, respectively. MCFs normally contain a central core that coincides with the fiber’s neutral axis, in addition to a number of symmetrically positioned outer cores. Gratings located within the central core are only sensitive to strain due to axial loading and temperature fluctuations; they are not sensitive to bending as the length of the neutral axis will not alter upon bending. Axial strain is generally assumed to be negligible, and thus the change in wavelength due to temperature change can be simply known by measuring the wavelength shift in the central core (). For the outer cores, the temperature contribution is thus known from the central core, and Eq. 1 can be rewritten as:where is the change in strain in the core. The relationship between the strain in each core and the corresponding curvature κ can be obtained from the geometry as depicted in Figure 1, and expressed as:where r is the distance of the outer cores to the central core assuming a symmetrical configuration (i.e. r is equal for all outer cores), and is the angle of the . core with respect to the x-axis, which intersects the second core. As can be seen from Eq. 3, the two unknowns are curvature κ and bend angle , where is defined as the angle between the x-axis and the curvature vector. Hence only two outer core strain measurements are needed to solve for them. However, additional cores can be used to improve the result and reduce errors. A closed-form solution can be obtained for the curvature κ and the bend angle by defining an apparent curvature vector ; :where and are the unit vectors along the x- and y-axes respectively, N is the number of outer cores, and represents the vector angle. The bend angle in Eq. 6 is obtained using the function. If a MCF contains n FBGs located at discrete intervals along its centerline, Eqs. 5 and 6 have to be computed n times. This results in a set of curvatures and bend angles along the fiber’s length with .
FIGURE 1
2.2 Cosserat Rod Model
The Cosserat rod theory presents a nonlinear mechanics-based model describing elastic rod deformations with regards to internal and external forces. Unlike simpler continuum robot modelling techniques, often relying on constant curvature assumptions, the Cosserat rod model provides accurate and geometrically exact solutions, even for large deflections and curvatures. Figure 2 shows an example of an elastic rod subjected to a variety of external forces. The model is characterized by a set of equilibrium equations and a set of constitutive equations.
FIGURE 2
2.2.1 Equilibrium Equations
These equations describe the relationship between external forces, namely distributed forces and distributed moments, and internal reactions, namely internal forces and internal moments. The equilibrium equations which are given as a set of first order differential equations and parametrized by the arc length s can be written as :where and are internal force and internal moment vectors, respectively, and are the angular and linear strain matrices in skew-symmetric form, respectively, and and are the external distributed force and distributed moment vectors, respectively. Note that all frames and vectors are given with respect to the rod’s material frame (also known as body frame) with their origins being along the rod’s centreline, see Figure 2.
2.2.2 Constitutive Equations
These equations describe the relationship between internal forces and moments with linear and angular strains as:where and are the linear and angular strain vectors, respectively, is the shear and axial stiffness matrix, is the bending and torsional stiffness matrix, and and are the reference linear and angular strain vectors, respectively. Note that the angular strain in Eq. 10 corresponds to:where κ is the curvature given in Eq. 5 and is the bend angle given in Eq. 6.
2.2.3 Boundary Conditions
The following assumptions made in this work are outlined first:
• the rod’s axial stiffness is much larger than the rod’s bending stiffness
• there is no, or negligible, linear strain, i.e. , such that axial forces are disregarded
• there is no, or negligible torsion, i.e.
• external loads are limited to point forces and distributed forces; there are no point moments nor distributed moments applied on the rod
The equilibrium and constitutive sets of equations are commonly solved using one of two constraint-based approaches: 1) an approach based on a boundary value problem (BVP), or 2) an approach based on an initial value problem (IVP). The BVP approach would require knowledge of the internal force and internal moment at the rod’s “boundaries”, i.e. at the base of the rod , and at the tip of the rod , where L is the rod’s length. However, the internal forces and moments are commonly unknown at the base. Iterative methods, e.g. shooting methods, are commonly employed to reduce the BVP into a system of solvable IVPs. Alternatively, considering that the rod is rigidly fixed, or cantilevered, at , it is known that the internal force and internal moment are both zero at if there is no external force/moment on the tip. Otherwise, the internal force/moment at the tip is equal to the external force/moment. Hence the problem can be reformulated as an IVP starting from with boundary conditions and (or and ), and solving backwards towards .
Notice that the set of equilibrium equations do not contain terms for externally applied point forces nor point moments . The way to include them is to divide the integration process into a smaller set of integrations based on the number and location of the external point loads, and including them as a boundary condition to the internal reactions at that location. For example, suppose that a rod is subjected to a point force and a point moment at locations and , respectively. The solution can then be found by dividing the problem into three integration sets, starting with and at , then at , and finally at . Here, the - and + indicates the value of or just before the given location (towards ) and just after the given location (towards ), respectively.
3 Body Force Estimation
In this section, two independent force estimation methods are proposed. The first group is based on piecewise polynomial segmentation of measured curvature, while the second group is based on model parameter estimation.
3.1 Piecewise Polynomial-Based Segmentation of Curvature
3.1.1 General Approach
Let us first consider the case where only point forces are applied onto the rod. Looking back at the Cosserat rod’s constitutive equations, it can be seen that the internal moment can be computed in Eq. 10 given knowledge of the angular strain vector . Notice that the stiffness matrix and reference angular strain are constants, and that can be obtained from the measured curvature as seen in Eq. 11. Accordingly, the internal moment computed in Eq. 10 can be numerically derived with respect to the arc length s. From the equilibrium equations, the derived internal moment can be replaced in Eq. 8 to compute the internal force (given that was assumed). Generally, the resulting force profile of should be continuous. However, when external loads with finite application length are applied, a discontinuity would occur. From this observation, the magnitudes and locations of the external point forces could be found at the break points (see Figure 3), i.e. points of discontinuity, in the internal force profile such that . The main drawback with this approach is that the curvature measurements obtained from the MCF tend to be somewhat noisy. Moreover, numerical derivation further amplifies the noise, which jeopardizes the quality of the force estimation.
FIGURE 3
3.1.2 Polynomials With n Degree
Let us now represent the curvature profile of a rod, subjected to a diversity of external forces, as a series of polynomials, each with its own degree n. According to the theory of elasticity, for a rod with constant flexural rigidity, the order of the curvature profile resulting from a given force type would be two degrees higher than the force profile itself. For example, a point force has a degree and would result in a linear curvature profile, . Similarly, a distributed force with constant magnitude has a degree and would result in a quadratic curvature profile . Hence, to overcome the noisy curvature measurements problem stated previously, the curvature profile can be segmented into a series of order polynomials for point loads, or a series of and order polynomials for distributed loads, or a combination of distributed loads and point loads. Note however, that for the final case, segmenting the curvature profile into both linear and quadratic segments is a very challenging task. Especially when there are sparse and discrete curvature measurements available, and when a quadratic segment can be very well approximated by a linear segment with negligible error. Accordingly, to prove the principle behind this approach, the estimation will be limited to point forces only, where an alternative technique will be provided in the second group of estimation methods for distributed forces.
3.1.3 Generic Concept
The concept behind piecewise polynomial-based segmentation of curvature, for the point forces case, is illustrated in Figure 3. Starting from a discrete set of noisy curvature measurements (green dots), a set of linear segments are to be constructed (blue lines) that resemble as close as possible the “ideal” theoretical curvature (dashed red lines). Several linear segmentation techniques can be found in literature, see for a review. The most common algorithms include: sliding window, top-down, and bottom-up. Qiao et al. previously employed a top-down linear segmentation approach to estimate one or two point forces applied on a Nitinol rod . The problem with top-down segmentation approaches, or most conventional linear segmentation approaches for that matter, is that they typically do not consider a penalty on the number of segments constructed. Furthermore, the top-down segmentation algorithm is known for its potential inflexibility in determining the break points .
3.2 Piecewise Polynomial-Based Segmentation Algorithms
Three alternative linear segmentation algorithms that are specific to the body force estimation problem and provide additional control where possible are proposed. All of the following algorithms start with raw curvature measurements that are merely interpolated to allow for a larger resolution along the arc length.
Algorithm 1
Curvature error thresholding
3.2.1 Curvature Error Thresholding
This curvature segmentation algorithm is illustrated in Algorithm 1, and the different properties are described in Table 1. The concept behind this method is that it iteratively varies the number of polynomial segments from , computes the error between the reconstructed segments and the measured curvature, and finally chooses the segments resulting in an error just below a predefined error .
TABLE 1
| Property | Definition |
|---|---|
| User defined curvature error threshold | |
| User defined maximum number of segments to iterate over | |
| Measured curvature from MCF | |
| Estimated force locations | |
| A sequence of interconnected linear curvature segments | |
| Structure to store the break points for each iteration; break points are locations of discontinuity in the curvature profile | |
| Structure to store the segmented curvature for each iteration | |
| Vector to store the error between and for each iteration | |
| User defined minimum detectable force | |
| User defined a priori estimate on the number of segments | |
| Minimum estimated force | |
| findBrkPnts (i, ) | Function to compute the optimal location of break points based on a desired number of break points i and a linear computational cost function, see for details |
| createSegments (, ) | Function to segment into a series of straight lines based on the given break points |
| indMinError (, ) | Function that starts from iteration , checks if the error in falls below ; if it does, it returns the number of break points according to the current iteration i, otherwise it increments i and repeats the process |
| findMinForce () | Function that propagates through the cosserat rod model and returns the minimum estimated force |
| mergeBrkPnts () | Function that merges the previously estimated break points based on the location where is found |
| mergeSegs () | Function that merges the previously estimated curvature segments based on the location where is found |
Definition of elements in the different algorithms.
3.2.2 Force Magnitude Thresholding
This curvature segmentation algorithm is illustrated in Algorithm 2, and the different properties are described in Table 1. The concept behind this method is that the measured curvature is initially segmented into a user-defined number of segments . These segments are then used to compute the forces using the Cosserat rod model. If a force magnitude is lower than a predefined minimum magnitude (this happens when the change in slope is very small), the corresponding two segments that caused this force are merged into one segment. The process repeats until there are no more forces below the predefined threshold.
3.2.3 Curvature Slope Change Thresholding
This method relies on the change point detection algorithm used in the findBrkPnts () function, see . This time however, instead of defining the number of segments to find the break points, a minimum slope change is provided. The algorithm consequently finds the number and location of break points, and the curvature is segmented accordingly.
For the three aforementioned methods, the output curvature segments are used in the constitutive Eq. 10 to compute the internal moment . The procedure then continues as explained previously in Section 3, A.1 to compute the magnitudes and locations of the applied external forces. Note that while the illustrated algorithms are restricted to linear curvature segmentation, they can be readily extended to combine both linear and quadratic curvature segmentation.
Algorithm 2
Force magnitude thresholding
3.3 Model-Based Estimation
Here, a predefined number N of forces and locations are estimated, where N must be larger or equal to the actual number of applied forces . If it is strictly larger, it will become clear that the magnitudes of those forces additional to will be equal to zero, or have a very small value. Let us now consider a function that requires an input and returns an output such that:where and are the magnitudes of the x and y components of the external point force relative to the rod’s body frame, respectively. is the location of the external point force expressed as arc length . and are the magnitudes of the x and y components of the curvature measurement along the arc length, and M the total number of curvature measurements over the total length. The function provides the forward kinematics of the Cosserat rod model, which is basically the inverse of the process described in Section 3, A. The estimation objective is to find given measured and forward kinematics function .
3.3.1 Least-Squares Optimization
An optimization algorithm typically finds a set of optimization parameters based on minimizing a desired cost function . In this case, the cost function is defined as:where is the measured curvature vector using the MCF. Hence, the goal is to minimize the error between the measured curvature and the estimated curvature , defined by the optimization parameters . While there is a myriad of optimization algorithms available in the literature, the interior, bounds-constrained trust region algorithm for nonlinear optimization is employed .
3.3.2 Unscented Kalman Filter (UKF)
In a similar way to optimization algorithms, Kalman filters can be used to estimate model parameters. In this case, the model parameters are the point force magnitudes and locations defined in . Compared to the linear Kalman filter (KF) or the Extended Kalman filter (EKF), the UKF is advantageous as it can accurately estimate the true mean and covariance of the states, in addition to the posterior mean and covariance to a third order for any nonlinearity. The EKF can only achieve first order accuracy. In addition, the computational complexity of the UKF is in the same order as that of the EKF . Finally, the UKF avoids the need to compute the state and measurement function jacobians, which sometimes have to be computed numerically.
The state transition model and the measurement model of the UKF are described as:where are the states at time step k such that , is the measurement function such that , is the measurement at time step k, and and are the process and measurement noise, respectively. and are the process and measurment noise covariance matrices at time step k, and are found by manual tuning. Note that the states , i.e. the magnitudes and locations of the point loads, are assumed to evolve through random walk. The remainder of the UKF equations follow the process outlined in , and are omitted here for brevity.
4 Experiments
4.1 Description of the Experimental Setup
Figure 4 illustrates the force estimation experimental setup. The setup comprises a Nitinol rod Figure 4 ⑥ with an outer diameter of 2.16 mm and an inner diameter of 1.65 mm. The rod is horizontally clamped to an aluminium frame using a chuck with a through hole. A low-friction PTFE tube with outer diameter of 1.60 mm and inner diameter of 0.25 mm is inserted into the lumen of the Nitinol rod. Accordingly, a MCF fiber with an outer diameter of 0.2 mm is inserted into the PTFE tube. At the proximal end, the PTFE tube is rigidly fixed to the Nitinol rod using adhesives. Similarly, the MCF is also rigidly fixed to the PTFE tube using adhesives. The Nitinol rod and PTFE tube combination have a flexural rigidity of 0.048 Nm2. The MCF (FBGS International NV, Geel, Belgium) has one central core, three outer cores, and a total of 18 FBG sets with a 10 mm FBG separation distance (see in Figure 1) covering a total length of 170 mm. Cable suspended weights Figure 4 ⑤ are used to apply point forces onto the rod in the direction of gravity. An ATI Nano17 six DoF force sensor (ATI Industrial Automation, North Carolina, United States) Figure 4 ③ is attached to a cable, which is wrapped around the Nitinol rod. This construction allows applying a known force (measured by the ATI nano17) in an arbitrary direction on the rod. An NDI electromagnetic tracking (EM) system (Northern Digital, Waterloo, Canada) Figure 4 ① is employed to determine the transformation between the cable tension force measured in the Nano17 sensor frame and the rod’s body frame. One EM tracking sensor Figure 4 ② is thus attached collinearly with the tensioned cable, and another is attached to the base of the rod. An optical fanout and interrogation system (FBGS International, Geel, Belgium) Figure 4 ⑦ is used to measure the wavelength shifts of the FBG sets within the MCF.
FIGURE 4
4.2 Calibration of EM Tracking Sensors With ATI Force Sensor
An EM tracking sensor is rigidly fixed to the cable attached to the ATI force sensor. The total applied force onto the Nitinol rod is equal to the tension force within this cable such that , where is the cable tension force as shown in Figure 4. Note that the weight of the cable is negligible compared to the applied forces. The z-axis of the EM tracking sensor frame is aligned with the cable and rigidly fixed in that configuration. To calibrate the direction of with respect to the rod’s body frame, two additional EM tracking sensors are used. The first is rigidly attached to the base of the rod at , and the second is attached at the tip of the rod at . The tip of the rod is then moved (by hand) in a spiral path and both EM tracking sensor data and the rod shape are recorded [see for the method used to reconstruct the rod’s shape]. The rigid transformation matrix EM between the EM frames and the rod’s shape is found using the rigid Iterative Closest Point (ICP) registration algorithm. Any EM frame can now be transformed to the base EM frame and consequently to the rod’s frame (hence the tip EM sensor is no longer needed and is removed). The cable tension force , which is aligned with the z axis of the cable EM frame is then transformed to the base EM frame and consequently to the rod’s body frame.
4.3 Description of Experiments
The following parameters were used in the experiments: , , N, , and m−1. Note that is a normalized value. Furthermore, , = diag [repmat (0.01 m−1, 1, 2M)], and = diag [repmat (0.1 N, 0.1 N, 0.1 m, 1, N)] where the function repmat (, , ) returns a matrix that is a repetition of the input matrix for rows and columns, and the function diag () returns a diagonal matrix based on the elements of .
Experiments were carried out for one, two and three point forces being applied onto the rod at a given time. Furthermore, for each load configuration, the location of the forces were changed/interchanged along the arc length and their magnitude was varied between 05.25 N. The ATI force sensor attached to the cable was also rotated slowly around the rod’s centreline to vary the applied two-dimensional forces quasi-statically. The mean, maximum, and standard deviation of the change in force magnitudes with respect to the ATI nano17 frame were 0.8 N/s, 16.9 N/s and 0.9 N/s respectively. The algorithms were running using c++ code on a Intel(R) Core(TM) i7-8850H, 2.60 GHz dual core processor on an Ubuntu 16.04 operating system.
5 Results
In total, there were eight single force, six dual force, and three triple force configurations. For each configuration the ATI force sensor was rotated along the rod’s centreline to vary the lateral x- and y-components of the force. This resulted in a total of 12,908 recorded forces applied onto the rod for all configurations. Forces with magnitude less than (defined earlier in Section 4, C) are discarded as they are considered to be negligible in magnitude and rod deformation. If a method estimates the number of forces to be less than the actual number of forces, then the error for the extra forces are the force magnitudes themselves and the error for their locations are taken as zero. Table 2 shows a summary of the results for the five force estimation methods (three piecewise polynomial-based curvature segmentation methods and two model-based estimation methods). Note that the errors are defined as the absolute value of the difference between ground truth and estimated forces/locations.
TABLE 2
| Method | |||||
|---|---|---|---|---|---|
| Result | 1.a | 1.b | 1.c | 2.a | 2.b |
| Computation time | 12 ms | 12 ms | 12 ms | 30 ms | 8 s |
| Correct no. of forces | 6,089 (47.2%) | 10,158 (78.7%) | 6,597 (51.1%) | 9,845 (76.3%) | 8,842 (68.5%) |
| [N] | 0.346 (6.6%) | 0.422 (8.0%) | 0.372 (7.1%) | 0.326 (6.2%) | 0.385 (7.3%) |
| [N] | 0.162 (3.1%) | 0.294 (5.6%) | 0.149 (2.8%) | 0.213 (4.1%) | 0.204 (3.9%) |
| [M] | 0.01160 (6.8%) | 0.00910 (5.4%) | 0.00973 (5.7%) | 0.00785 (4.6%) | 0.00714 (4.2%) |
| [N] | 0.315 (6.0%) | 0.345 (6.6%) | 0.310 (5.9%) | 0.309 (5.9%) | 0.314 (6.0%) |
| [N] | 0.242 (4.6%) | 0.392 (7.5%) | 0.253 (4.8%) | 0.218 (4.2%) | 0.222 (4.2%) |
| [M] | 0.01551 (9.1%) | 0.00681 (4.0%) | 0.01040 (6.1%) | 0.00567 (3.3%) | 0.00588 (3.5%) |
Summary of force estimation results illustrating estimated force numbers and errors on magnitudes and locations. and represent the mean and standard deviation of the absolute value of the error, respectively.
It can be generally seen that the model-based estimation methods provide an overall better performance compared to their polynomial-based segmentation counterparts. This is because they can more often correctly estimate the actual number of forces, in addition to having a lower mean error on the force locations. Note however, that the force magnitude errors are comparable for all methods. The computational times of the polynomial-based segmentation methods are significantly lower however, taking around 1–2 ms for completion. The UKF approach takes around 30 ms, while the least-squares optimization method takes over 8 s. This is the major drawback of the least-squares optimization method, making it unusable for practical real-time settings.
Figure 5 shows an example result for a configuration where three external point forces are applied onto the rod. The forces at the extremities are due to suspended weights, while the force in the middle is due to the ATI force sensor and the attached cable. The curvature profiles in the x- and y-directions clearly illustrate how the different point forces in given directions correspond to changes of slope. It can also be seen that the measured curvature, which is a linear interpolation of the discrete curvature measurements at each FBG location, deviates from the ground truth curvature. This can be due to a variety of reasons including: having a gap between the MCF and the inner lumen of the surrounding tube, friction between the MCF and the inner lumen of the surrounding tube which causes twist during bending, and FBG sensor accuracy itself. These directly impact the measured curvature κ and bend angle profiles, which are the key to accurate force estimation. Figure 5 also shows the deformed Nitinol rod and a comparison between the ground truth shape, estimated shape, and the applied forces. While the measured curvature may be noisy for direct force estimation, it can be readily employed to reconstruct the shape of the rod, with tip errors 2 mm. Furthermore, the close correspondence between the ground truth forces and locations with the estimated forces and locations can also be clearly seen. Notice that in this configuration, an extra force is estimated (see Figure 5). However, because its magnitude is very small compared to the other forces and less than , it is assumed negligible and discarded. This is often the case since the model assumes four applied forces a priori, while in reality, only one to three forces are applied.
FIGURE 5
6 Distributed Forces
Thus far, the methods and experimental validations described concerned point force estimation. However, in some applications forces can be distributed. This section introduces some first works explaining how distributed forces can be estimated. Note that these methods concern distributed forces with constant magnitude.
6.1 Piecewise Polynomial-Based Segmentation of Curvature
It was previously pointed out that it is possible to estimate distributed forces using piecewise polynomial-based curvature segmentation. The method must be extended to allow for 1) second order polynomial detection, and 2) classification of polynomial segments to first or second degree. However, when curvature measurements are discrete, sparse, and require interpolation to provide extra smoothness (as is the case here since curvature measurements are only available at 10 mm intervals), it becomes very challenging to distinguish between first order and second order polynomials, especially if the curvature measurements are noisy as well. This task, however, can be simplified if extra information was given a priori such as the locations of the forces or the total number of applied forces. Since this information is not assumed available, we will not continue with this method for distributed forces estimation.
6.2 Model-Based Estimation
6.2.1 Method Description
The model-based method is more general and flexible as the underlying Cosserat rod theory is general and can be used to model the rod’s forward kinematics using any type of external force. Hence, the strategy to estimate distributed forces is a simple expansion to the method developed for point forces in Section 3B. The difference is in the model’s input vector , which is now defined as:where and are the magnitudes of the and components of the external distributed force, respectively, and and are the start and end locations of the distributed force along the arc length. The output curvature vector and the forward kinematics function remain the same. The only difference is that distributed forces were set to zero in the earlier implementation, whereas now they are not.
6.2.2 UKF Simulation
To illustrate the feasibility of the proposed method, a simulation study is provided using the UKF approach. Note that the least-squares optimization approach can be employed in theory, but is disregarded in this study for its high computational cost. Here, two distributed forces with varying magnitudes and locations are applied onto the rod. The magnitudes and locations vary according to a sinusoidal profile as shown in Figure 6A. Gaussian noise with zero mean and 0.01 m−1 standard deviation was added to the curvature measurements. The UKF estimation of the distributed force magnitudes and locations is also shown for a simulation time of 5 s. Table 3 gives a summary on the absolute value of the errors between the ground truth and estimated forces and locations. It is clear that the UKF shows excellent performance in estimating both their magnitudes and locations. When comparing the mean magnitude and location errors with the experimental errors shown in Table 2, it is evident they are lower, and that the overall performance is better. Even though the results shown are for a simulation study only, they do give an indication that the performance in practice could be decent as well. Finally, a similar simulation was carried out by applying one point force and one distributed force. The same distributed force UKF model was able to estimate the point force magnitude and location to the same accuracy as outlined previously. This shows that the distributed force model is a generic model allowing for the estimation of generic force types.
FIGURE 6
TABLE 3
| [N/m] | [N/m] | [m] | [m] | |
|---|---|---|---|---|
| Mean | 1.19 (1.6%) | 0.47 (0.6%) | 0.0018 (1.1%) | 0.0017 (1.0%) |
| Max | 10.21 (13.6%) | 13.73 (18.3%) | 0.0107 (6.3%) | 0.0147 (8.6%) |
Simulation distributed force estimation errors over time.
7 Conclusion and Discussion
Estimation of forces along instrument bodies is a crucial aspect in medical navigation and handling. However, estimating such forces has proven to be challenging. Sensitive tactile sensors could solve this problem, but designing them to be robust and compact, featuring high spatial and force resolution has been found extremely difficult. The recent advent of compact MCFs opens up some opportunities for distributed contact force sensing, as demonstrated in this work. Here the instrument’s deformation itself can be exploited to provide estimates of the external forces. Two approaches have been developed to provide point and distributed force estimations using curvature information exclusively from a MCF. The resulting performance, for all methods, clearly indicates their feasibility for body force estimation. The different methods offer different advantages and result in an average of around 64.4% correct number of forces estimations, 5.5% force magnitude error, and 5.3% location error.
In the latter methods, the model-based UKF has shown to stand out when comparing accuracy, generalizability and computational time of the different methods. Furthemore, it must be noted that the Nitinol rod used in the validation experiments has a relatively high flexural rigidity, EI = 0.048 Nm2. In reality, medical instruments tend to have lower flexural rigidities. This is more beneficial as it allows for larger curvatures and an easier distinction between consecutive segments. It also overcomes minimum sensitivity issues regarding the MCF. Hence the proposed methods will work better for instruments with relatively low flexural rigidities. However, a key aspect to keep in mind is the elasticity of the instrument, which is necessary to ensure repeatability of the results.
Finally, a generalized method to estimate externally applied distributed forces was proposed, and illustrated that it can actually be applied for both point forces and distributed forces. Extensions to the proposed model can be easily implemented for diverse types of external loads. Furthermore, the feasibility of the proposed method was proved for forces with varying magnitudes and locations. Meaning that the method not only works for static loads, but also for quasi-static loads. Future work will focus on validating this method experimentally on elastic instruments with low flexural rigidities.
Statements
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author contributions
In this work, OA and EVP proposed the research topic as an extension to the state-of-the-art methods for estimation of external forces. OA conducted the simulations, experiments, data post-processing and analysis, and writing of the manuscript. MO, JV, and EVP provided continuous guidance and feedback regarding the proposed force estimation strategies and brought about different ideas for the direction of the research. Furthermore, MO, JV, and EVP provided valuable feedback and review of the manuscript.
Funding
The work was funded by the Flemish Agency of Innovation and Entrepreneurship (VLAIO), under project number (HBC.2018.2046).
Conflict of interest
Authors OA and JV are employed by FBGS International NV.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
Al-AhmadO.OurakM.Van RoosbroeckJ.VlekkenJ.Vander PoortenE. B. (2020). Improved FBG-Based Shape Sensing Methods for Vascular Catheterization Treatment. IEEE Robot. Autom. Lett.5, 1. 10.1109/LRA.2020.3003291
2
AloiV. A.RuckerD. C. (2019). “Estimating Loads along Elastic Rods,” in IEEE International Conference on Robotics and Automation, Montreal, QC, Canada, May 20–24, 2019, 2867–2873. 10.1109/ICRA.2019.8794301
3
AntmanS. S. (2005). Nonlinear Problems of Elasticity. Appl. Math. Sci.107, 1–835. 10.1007/978-1-4757-4147-6_1
4
BackJ.ManwellT.KarimR.RhodeK.AlthoeferK.LiuH. (2015). “Catheter Contact Force Estimation from Shape Detection Using a Real-Time Cosserat Rod Model,” in 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Hamburg, Germany, September 28–October 2, 2015, 2037–2042. 10.1109/IROS.2015.7353647
5
BalamiJ. S.WhiteP. M.McMeekinP. J.FordG. A.BuchanA. M. (2018). Complications of Endovascular Treatment for Acute Ischemic Stroke: Prevention and Management. Int. J. Stroke13, 348–361. 10.1177/1747493017743051
6
Burgner-KahrsJ.RuckerD. C.ChosetH. (2015). Continuum Robots for Medical Applications: A Survey. IEEE Trans. Robot.31, 1261–1280. 10.1109/TRO.2015.2489500
7
ColomeA.PardoD.AlenyaG.TorrasC. (2013). “External Force Estimation during Compliant Robot Manipulation,” in 2013 IEEE International Conference on Robotics and Automation, Karlsruhe, Germany, May 6–10, 2013, 3535–3540. 10.1109/ICRA.2013.6631072
8
DalyJ. M.WangD. W. L. (2014). Time-Delayed Output Feedback Bilateral Teleoperation with Force Estimation for $n$-DOF Nonlinear Manipulators. IEEE Trans. Contr. Syst. Technol.22, 299–306. 10.1109/TCST.2013.2242329
9
De GreefA.LambertP.DelchambreA. (2009). Towards Flexible Medical Instruments: Review of Flexible Fluidic Actuators. Precision Eng.33, 311–321. 10.1016/j.precisioneng.2008.10.004
10
GallowayK. C.ChenY.TempletonE.RifeB.GodageI. S.BarthE. J. (2019). Fiber Optic Shape Sensing for Soft Robotics. Soft Robot.6, 671–684. 10.1089/soro.2018.0131
11
GaoA.LiuN.ShenM.E.M.K. AbdelazizM.TemelkuranB.YangG.-Z. (2020). Laser-Profiled Continuum Robot with Integrated Tension Sensing for Simultaneous Shape and Tip Force Estimation. Soft Robot.7, 421–443. 10.1089/soro.2019.0051
12
HaideggerT.BenyóB.KovácsL.BenyóZ. (2009). Force Sensing and Force Control for Surgical Robots. IFAC Proc. Volumes42, 401–406. 10.3182/20090812-3-DK-2006.0035
13
HasanzadehS.Janabi-SharifiF. (2015). Model-based Force Estimation for Intra-cardiac Catheters. IEEE ASME Trans. Mechatron.21, 1. 10.1109/TMECH.2015.2453122
14
HeunisC. M.BelfioreV.VendittelliM.MisraS. (2019). “Reconstructing Endovascular Catheter Interaction Forces in 3D Using Multicore Optical Shape Sensors,” in In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Macau, China, November 3–8, 2019, 5419–5425. 10.1109/IROS40897.2019.8968448
15
HeunisC.SikorskiJ.MisraS. (2018). Flexible Instruments for Endovascular Interventions: Improved Magnetic Steering, Actuation, and Image-Guided Surgical Instruments. IEEE Robot. Automat. Mag.25, 71–82. 10.1109/MRA.2017.2787784
16
HooshiarA.SayadiA.JolaeiM.DargahiJ. (2020). “Accurate Estimation of Tip Force on Tendon-Driven Catheters Using Inverse Cosserat Rod Model*,” in 2020 International Conference on Biomedical Innovations and Applications (BIA), Varna, Bulgaria, September 24–27, 2020, 37–40. 10.1109/BIA50171.2020.9244512
17
HuJ.XiongR. (2018). Contact Force Estimation for Robot Manipulator Using Semiparametric Model and Disturbance Kalman Filter. IEEE Trans. Ind. Electron.65, 3365–3375. 10.1109/TIE.2017.2748056
18
JungJ.LeeJ.HuhK. (2006). Robust Contact Force Estimation for Robot Manipulators in Three-Dimensional Space. Proc. Inst. Mech. Eng. C. J. Mech. Eng. Sci.220, 1317–1327. 10.1243/09544062C09005
19
KavicS. M.BassonM. D. (2001). Complications of Endoscopy. Am. J. Surg.181, 319–332. 10.1016/S0002-9610(01)00589-X
20
KhanF.RoesthuisR. J.MisraS. (2017). “Force Sensing in Continuum Manipulators Using Fiber Bragg Grating Sensors,” in 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Vancouver, BC, Canada, September 24–28, 2017, 2531–2536. 10.1109/IROS.2017.8206073
21
KhoshnamM.SkanesA. C.PatelR. V. (2015). Modeling and Estimation of Tip Contact Force for Steerable Ablation Catheters. IEEE Trans. Biomed. Eng.62, 1404–1415. 10.1109/TBME.2015.2389615
22
KillickR.FearnheadP.EckleyI. A. (2012). Optimal Detection of Changepoints with a Linear Computational Cost. J. Am. Stat. Assoc.107, 1590–1598. 10.1080/01621459.2012.737745
23
LeT. M.FatahiB.KhabbazH.SunW. (2017). Numerical Optimization Applying Trust-Region Reflective Least Squares Algorithm with Constraints to Optimize the Non-linear Creep Parameters of Soft Soil. Appl. Math. Model.41, 236–256. 10.1016/j.apm.2016.08.034
24
LovrićM.MilanovićM.StamenkovićM. (2014). Algoritmic Methods for Segmentation of Time Series: An Overview. J. Contemp. Econ. Bus. Issues1, 31–53.
25
LvX.WuZ.JiangC.LiY.YangX.ZhangY. (2011). Complication Risk of Endovascular Embolization for Cerebral Arteriovenous Malformation. Eur. J. Radiol.80, 776–779. 10.1016/j.ejrad.2010.09.024
26
MassariL.SchenaE.MassaroniC.SaccomandiP.MenciassiA.SinibaldiE.et al (2020). A Machine-Learning-Based Approach to Solve Both Contact Location and Force in Soft Material Tactile Sensors. Soft Robot.7, 409–420. 10.1089/soro.2018.0172
27
MooreJ. P.RoggeM. D. (2012). Shape Sensing Using Multi-Core Fiber Optic cable and Parametric Curve Solutions. Opt. Express20, 2967–2973. 10.1364/oe.20.002967
28
QiaoQ.BorghesanG.De SchutterJ.Vander PoortenE. (2021). Force from Shape [start]â[end][start]€[end][start]”[end] Estimating the Location and Magnitude of the External Force on Flexible Instruments. IEEE Trans. Robot.
29
QiaoQ.WillemsD.BorghesanG.OurakM.De SchutterJ.PoortenE. V. (2019). “Estimating and Localizing External Forces Applied on Flexible Instruments by Shape Sensing,” in 2019 19th International Conference on Advanced Robotics (ICAR), Belo Horizonte, Brazil, December 2–6, 2019, 227–233. 10.1109/ICAR46387.2019.8981592
30
RuckerD. C.WebsterR. J. (2011). “Deflection-based Force Sensing for Continuum Robots: A Probabilistic Approach,” in 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Francisco, CA, September 25–30, 2011, 3764–3769. 10.1109/IROS.2011.6048202
31
RyuH. T.WooJ.SoB. R.YiB. J. (2020). Shape and Contact Force Estimation of Inserted Flexible Medical Device. Int. J. Control Autom. Syst.18, 163–174. 10.1007/s12555-019-0237-8
32
SadatiS. M.ShivaA.HerzigN.RuckerC. D.HauserH.WalkerI. D.et al (2020). Stiffness Imaging with a Continuum Appendage: Real-Time Shape and Tip Force Estimation from Base Load Readings. IEEE Robot. Autom. Lett.5, 2824–2831. 10.1109/LRA.2020.2972790
33
ShiC.LuoX.QiP.LiT.SongS.NajdovskiZ.et al (2017). Shape Sensing Techniques for Continuum Robots in Minimally Invasive Surgery: A Survey. IEEE Trans. Biomed. Eng.64, 1665–1678. 10.1109/TBME.2016.2622361
34
TrejosA. L.PatelR. V.NaishM. D. (2010). Force Sensing and its Application in Minimally Invasive Surgery and Therapy: A Survey. Proc. Inst. Mech. Eng. Part. C J. Mech. Eng. Sci.224, 1435–1454. 10.1243/09544062JMES1917
35
VenkiteswaranV. K.SikorskiJ.MisraS. (2019). Shape and Contact Force Estimation of Continuum Manipulators Using Pseudo Rigid Body Models. Mech. Mach. Theor.139, 34–45. 10.1016/j.mechmachtheory.2019.04.008
36
WanE. A.Van Der MerweR. (2000). “The Unscented Kalman Filter for Nonlinear Estimation,” in Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium (Cat. No.00EX373), Lake Louise, AB, Canada, October 4, 2000, 153–158. 10.1109/ASSPCC.2000.882463
37
WuJ.LiB.LiC.QiuY.ZhangZ.LyuW.et al (2021). “Fiber Bragg Grating-Based Shape Sensing: a Review and Perspective,” in Global Intelligent Industry Conference 2020, Guangzhou, China, November 16–18, 2020. 10.1117/12.2590684
38
XuK.SimaanN. (2010). Intrinsic Wrench Estimation and its Performance index for Multisegment Continuum Robots. IEEE Trans. Robot.26, 555–561. 10.1109/TRO.2010.2046924
39
YasinR.SimaanN. (2020). Joint-level Force Sensing for Indirect Hybrid Force/position Control of Continuum Robots with Friction. Int. J. Rob. Res.10.1177/0278364920979721
Summary
Keywords
external force estimation, fiber bragg gratings, multi-core fibers, cosserat rod, flexible instruments
Citation
Al-Ahmad O, Ourak M, Vlekken J and Vander Poorten E (2021) FBG-Based Estimation of External Forces Along Flexible Instrument Bodies. Front. Robot. AI 8:718033. doi: 10.3389/frobt.2021.718033
Received
31 May 2021
Accepted
09 July 2021
Published
30 July 2021
Volume
8 - 2021
Edited by
Chaoyang Shi, Tianjin University, China
Updates
Copyright
© 2021 Al-Ahmad, Ourak, Vlekken and Vander Poorten.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Omar Al-Ahmad, omar.alahmad@kuleuven.be
This article was submitted to Biomedical Robotics, a section of the journal Frontiers in Robotics and AI
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.