A Three-Dimensional Shape-Based Force and Stiffness-Sensing Platform for Tendon-Driven Catheters

This paper presents an efficient shape-based three-axial force and stiffness estimator for active catheters commonly implemented in cardiac ablation. The force-sensing capability provides important feedback for catheterization procedures including real-time control and catheter steering in autonomous navigation systems. The proposed platform is based on the introduced accurate and computationally efficient Cosserat rod model for tendon-driven catheters. The proposed nonlinear Kalman filter formulation for contact force estimation along with the developed catheter model provides a real-time force observer robust to nonlinearities and noise covariance uncertainties. Furthermore, the proposed platform enables stiffness estimation in addition to tip contact force sensing in different operational circumstances. The approach incorporates pose measurements which can be achieved using currently developed pose-sensing systems or imaging techniques. The method makes the approach compatible with the range of forces applied in clinical applications. The simulation and experimental results verify the viability of the introduced force and stiffness-sensing technique.


Introduction
Many medical interventions such as cardiac interventions and interventional radiology involve deployment of steerable catheters for diagnosis and treatment of internal body organs [1,2]. Despite recent advances in robot-assisted catheterization [2][3][4][5][6], the catheterization is usually performed manually by rotating and translating the proximal end of the catheter as well as bending the distal shaft via a steering knob on the handle (Figure 1). Steering of catheters, due to their inherent flexibility and operation in a confined environment, requires high dexterity and expertise to navigate and place the tip on the target tissue. Tip force measurement could have great implications in the success of both manual and robot-assisted procedures. For instance, during the cardiac ablation, a compliance force should be applied by the end point to the intra-cardiac tissue for a certain period of time and the tip of catheter must maintain a stable contact with the target, whereas inadequate force monitoring results in complications, including perforation and permanent destruction of cardiac tissues due to excessive force application or recurrence of the cardiac problem owing to low force implementation [1,7,8].
Despite its importance, force feedback remains one of the foremost challenges in advancing catheterization techniques. To this end, two commercial tri-axial force-sensing ablation catheters have been recently introduced, namely TactiCath ® catheter system (St. Jude Medical, St. Paul, MN, USA) and Thermocool ® SmartTouch™ (Biosense Webster, Diamond Bar, CA, USA) [9,10]. The TactiCath ® have been recently introduced, namely TactiCath ® catheter system (St. Jude Medical, St. Paul, MN, USA) and Thermocool ® SmartTouch™ (Biosense Webster, Diamond Bar, CA, USA) [9,10]. The TactiCath ® is a 3.5-mm open-irrigated tip radiofrequency (RF) ablation catheter which uses tri-axial fiberoptic cable for contact force sensing. The Thermocool ® SmartTouch™ has a 3.5-mm six-hole open-irrigated tip that measures the force by deformation of a precision spring which is located between the tip and distal shaft of the catheter. Both catheters measure and visualize the magnitude and direction of tip contact forces, but they require accompanying systems to operate with. Therefore, equipping the operating room with the specialized navigation system in addition to the higher cost of the force-sensing catheters impose significant limitations to their use. Some recent studies have focused on the development of sensors and catheters with force sensors, including Polygerinos et al., who developed three-dimensional magnetic resonance imaging (MRI)-compatible fiber optic force sensor with the reported force resolution of 0.01 N and 5% error in force measurement with the maximum magnitude of 0.5 N [11]. Kenser and Howe developed an axial fiber optic force sensor and mounted it on the ablation catheter tip [12,13]. The prototype is approximately 15Fr (5 mm) in diameter with 2%-4% error in axial force sensing.
A multi-axis force sensor at the tip of the tool directly measures the interaction force. However, the introduced sensors suffer from significant impediments. They are limited in the number of measurement axes, and introduce extra cost, packaging issues, sensitivity to temperature, and limitations in MRI compatibility. Integration of the sensor with conventional systems may also lead to interference with the operation of the ablation catheters, and create an extra obstacle to the effort for producing smaller end-effectors. To cope with these problems, the tip-tissue contact force would ideally be estimated without using force sensors.  Xu et al. presented the idea of "intrinsic force sensing" for continuum robots that is based on the evaluation of certain components of tip forces through the joint actuation loads [14], and introduced performance index for analysis of full-wrench sensing [15]. Inspired by the suggested concept, Rucker and Webster proposed a probabilistic approach for planar continuum robot force estimation using tip pose and actuation torques [16]. However, their method is limited to planar forces and suffers from convergence issues. Khoshnam et al. investigated contact force estimation with the Euler-Bernouli model for the catheter deflection [17], where the planar contact forces were estimated using a strain sensor mounted at the distal base and numerical mapping of the tip angles to the sensor readings. In a recent study [18], it was demonstrated that the flexibility of a catheter's distal shaft can be implemented as a means of tip contact force estimation. To this end, the same group introduced the force index and provided only the contact force range with 80% accuracy. In summary, the aforementioned approaches have one or more of the following limitations: assumption of planar forces/deformation, small deformations, and constant curvature or piece-wise constant curvature, in addition to convergence issues, sensitivity to noise, and need for actuation torque measurements.
To address the aforementioned problems, this paper contributes by developing a 3D force estimation platform for tendon-driven catheters. For this purpose, we introduced the more powerful Xu et al. presented the idea of "intrinsic force sensing" for continuum robots that is based on the evaluation of certain components of tip forces through the joint actuation loads [14], and introduced performance index for analysis of full-wrench sensing [15]. Inspired by the suggested concept, Rucker and Webster proposed a probabilistic approach for planar continuum robot force estimation using tip pose and actuation torques [16]. However, their method is limited to planar forces and suffers from convergence issues. Khoshnam et al. investigated contact force estimation with the Euler-Bernouli model for the catheter deflection [17], where the planar contact forces were estimated using a strain sensor mounted at the distal base and numerical mapping of the tip angles to the sensor readings. In a recent study [18], it was demonstrated that the flexibility of a catheter's distal shaft can be implemented as a means of tip contact force estimation. To this end, the same group introduced the force index and provided only the contact force range with 80% accuracy. In summary, the aforementioned approaches have one or more of the following limitations: assumption of planar forces/deformation, small deformations, and constant curvature or piece-wise constant curvature, in addition to convergence issues, sensitivity to noise, and need for actuation torque measurements.
To address the aforementioned problems, this paper contributes by developing a 3D force estimation platform for tendon-driven catheters. For this purpose, we introduced the more powerful Cosserat rod-based quasi-static model of tendon-driven catheters that can be used for both small and large 3D deformations. The computational efficiency of the model along with its accuracy makes it appealing for real-time applications. Based on the model, tri-axial shape-based contact force and stiffness estimators for tendon-driven catheters are developed which do not require actuation loading information. Potential for stiffness estimation is another advantage of the model. Tip pose and shape measurements, readily available in clinical applications through imaging, can be deployed for tip contact force estimation avoiding extra sensor integration [19,20]. In this study, electromagnetic tracking sensors were incorporated for pose measurements. Several force estimation platforms based on the Kalman filter (KF) are designed and their performances are compared in terms of efficiency, convergence, and robustness. The adoption of a KF-based approach makes it less sensitive to noise and enables next time step force prediction (useful for tracking control purposes).

Modeling the Steerable Catheter
A model of force displacement is required for shape-based force estimation of tendon-driven catheters. Therefore, with the introduction of the Cosserat theory of elastic rods [21][22][23], we develop the quasi-static model for tendon-driven intra-cardiac catheters. For convenience, the nomenclature is introduced in Table 1.

Kinematics
An elastic rod with an unstretched length L is characterized by the global coordinate frame e " te 1 , e 2 , e 3 u represented in Figure 2a. The vector r(s), s P r0, Ls locates the centerline of the rod with respect to the global coordinate. The local frame dpsq " td 1 psq, d 2 psq, d 3 psqu attached to the curve along the centerline demonstrates the orientation R c psq P SOp3q of each point on the rod relative to the global frame. The director d 3 is adopted to be tangent to the centerline at each infinitesimal segment of the rod and d 1 and d 2 define the local frame on the cross-section of the rod. The velocity vector v l " rv l 1 , v l 2 , v l 3 s T defines shear strain along d 1 and d 2 by the first two components, and axial compression or stretch along d 3 by v l 3 , where superscript l indicates definition with respect to the local frame. With the shear strain and elongation being negligible, the change of the rod's position with respect to arc length becomes, The rate of change of each director along arc length is, where the Darboux vector is given by ω " rk 1 , k 2 , τs T with k i pi " 1, 2q, denoting curvature along d i and τ representing the torsion.
and representing the torsion.

Governing Equilibrium Equations
The governing equilibrium equation of model in global frame is described by where dot represents the derivative with respect to arc length [22]. The internal moment and force vectors, represented by l l l lT s n n n  respectively, are described with respect to the local frame. Internal moments are related to the local kinematic displacement of the rod as follows,

Applied External Loads to the Ablation Catheter
The primary interacting forces for a tendon-driven catheter include tendon tension and cardiac-interacting load.

Governing Equilibrium Equations
The governing equilibrium equation of model in global frame is described by where dot represents the derivative with respect to arc length [22].
The internal moment and force vectors, represented by m l psq " rm l 1 m l 2 m l 3 s T and n l psq " rn l 1 n l 2 n l 3 s T , respectively, are described with respect to the local frame. Internal moments are related to the local kinematic displacement of the rod as follows, m l 1 psq " EI 1 pk 1 psq´k 1,0 q, m l 2 psq " EI 2 pk 2 psq´k 2,0 q, m l 3 psq " GJ s pτpsq´τ 0 q (4)

Applied External Loads to the Ablation Catheter
The primary interacting forces for a tendon-driven catheter include tendon tension and cardiac-interacting load.
The tendon tension force, F T , leads to external torque T l " F T d tc at the tip along d 1 resulting from the distance d tc between the tendon and centerline (Figure 1b). This initial torque can be modeled as an initial curvature, k 1,0 , as well [24,25].
During the ablation, consistent contact with the cardiac wall is required. The contact changes as a result of heartbeat and respiration that leads to deviation in the contact forces. Assuming that the contact force is not available, the shape and deflection of the catheter should be measured continuously to give the estimate of changing forces.
The underlying assumptions are as follows, 1. Gravity effect due to the low mass and inertia of distal shaft is ignored. 2. Blood drag force is negligible when compared to the tip contact and tension forces [23,26]. 3. Tendon actuation effect is modeled as an end-moment load implemented at the tip of the distal shaft. In this regard, it is assumed that tendon's axial tip force and distributed forces cancel each other, and the range of transverse tip forces and curvatures do not lead to significant model inaccuracy [23,27,28].
Therefore, we assume fpsq " lpsq " 0, and consider the external point forces F to be applied at the tip of the rod.

Quasi-Static Model of Tendon-Driven Catheter
By substituting Equation (2) into Equation (3), the resulting internal force equilibrium equation is obtained. For ease of notation, dependence on s is dropped hereafter. Therefore, we have . n l 1`k 2 n l 3´τ n l 2 " 0 . n l 2`τ n l 1´k 1 n l 3 " 0 . n l 3`k 1 n l 2´k 2 n l 1 " 0 Likewise, replacing Equation (2) into Equation (3), and deploying Equations (1) and (4) leads to: The contact forces n l plq " rF x F y F z s T and contact torques m l plq " rT 1 T 2 T 3 s T are the initial boundary values of the above differential equations.
The rotation matrix can be specified using Euler angle representation by three successive rotations; first, around the e 2 axis by an angle of ψ, then the currentê 1 axis by θ, and finally the d 3 axis by ϕ. Therefore, » -- Numerical integration of Equations (7) and (8) specifies orientation and spatial position of all points and particularly the end point, required in the following section, on the centerline of the catheter. Equations (5) and (6) as an Initial Value Problem (IVP) provide pk 1 , k 2 , τq to be substituted in Equation (7).

Kalman Filter-Based Force Observer Design
In many interventions, there are noise contaminations. Moreover, system dynamics and measurement random process properties may vary. Real-time estimation of force imposes another constraint on the efficiency of the solution. Therefore, analytical solutions might not necessarily provide effective solutions. Significantly, stiffness of the catheter might not be available in advance or it could even change during the operation. Kalman filter-based estimators do not require deterministic dynamics or accurate measurements or the random process with stationary properties. In addition to force, stiffness can also be included in the state vector for real-time estimation. Additionally, the noise-filtering advantage of Kalman filters (KFs), as well as their inherent simplicity and effectiveness make KFs well accredited for real applications [29].
In this section, we develop a force estimator based on the adapted model, and compare variants of introduced KF-based estimators with respect to their performance. The proposed scheme allows us not only to provide contact force estimations but also to accommodate the estimates of dominant features causing a catheter's deformation, e.g., stiffness.

Contact Force Estimation (CFE)
For catheter contact force sensing, the desired variables are three-dimensional forces exerted on the catheter tip in local frame. Therefore, we define the state vector of the estimator as follows, It has been shown that contact forces can be observed from the catheter shape and deflection [18]. The curvatures, torsion and forces of the segments of the catheter are related through Equations (5) and (6). Consequently, according to Equations (7) and (8), catheter tip pose at time k, TpX k q " rx e,k y e,k z e,k θ e,k ψ e,k ϕ e,k s T , is a function of applied tip forces. Additionally, regarding corresponding equations, it is concluded that contact forces alter the local curvatures and torsion. Hence, the state space model of the estimator can be formulated as follows, where the catheter's current state, X k , is assumed to be related to the previous state, X k´1 , with model uncertainties. The sensor measurements, Y k , comprises tip pose and the second-last segment deformations ω se " rk 1,se k 2,se τ se s T as indicated in Figure 2. Here, η k and w k are independent zero-mean Gaussian noises that denote uncertainty in states and measurements. Also hpX k q " rTpX k q ω se pX k qs T By exploiting the developed model, the Jacobian matrix J " rJ TF J ωF s T is obtained through calculation of pose (T) and deformation (ω se ) changes with respect to contact force variations. These matrices can be computed both numerically through approximation with finite differences or analytically with higher accuracy at the expense of more complexity [30].

Contact Force and Stiffness Estimation (CFSE)
In order to estimate deflectable distal shaft's stiffness along with contact force, the contact force and bending stiffness are assumed as the state vector of the estimator as follows: where v p is the Poisson's ratio. The measurements are considered the same as those presented in Equation (11). Hence, incorporating the state space model specified in Equation (10), the desired states can be estimated. The Jacobian matrix J a " rJ Ta J ωa s T is determined through calculation of pose (T) and deformation (ω se ) changes with respect to contact force and stiffness variations. Therefore, J Ta " are established by the definition of J TE and J ωE as:

Extended Kalman Filtering
In practice, to deal with system nonlinearities (of state transition and observation functions), extended KF (EKF) is used in many engineering applications [29]. The filter is based on linearization of the nonlinear model with the first order Taylor expansion.
The EKF-based estimator for both CFE and CFSE is formulated as follows, HereX k , P k and K k denote the estimated state for both CFE and CFSE, error covariance matrix and Kalman gain at each time step, respectively. The initial values of the system state vector and error covariance matrix (measured through experiments or initial guess), are represented byX 0 and P 0 , respectively. Significant errors in the initial state estimate, process model, linearization, and covariance matrices may lead to divergence or instability of the estimates in EKFs.

Iterative and Adaptive Kalman Filters
In order to cope with linearization errors, iterative Kalman filters have been proposed [31], while adaptive filters have been introduced to deal with noise covariance matrix uncertainties [32]. These approaches are concisely presented which can be hierarchically deployed into the structure of nonlinear KF-based estimators.

Iterative Schemes
The iterative filter uses measurements in a recurrent loop to approximate the nonlinear function closely. Additionally, the scheme mitigates the sensitivity to error covariance matrix [31]. More precisely, the prediction equations of the EKF-based estimators are the same in this scheme.
However, the a posteriori state estimateX where l i is the number of iteration, is repeatedly calculated using the state estimate of previous iterations. At the end of iterationsX k "X N k , and the a posteriori covariance matrix estimate is updated withX k .

Adaptive Schemes
To deal with uncertainties of model and measurements, adaptive filters are introduced [32][33][34]. Various types of adaptive filters can be adopted. Nonetheless, since the results of different variants are similar, covariance matching is adopted. The algorithm is based on computing noise statistics for the estimator. The algorithm starts with an initial guess forq whereq k andr k are the state and measurement noise samples, respectively. These estimates with the noise covariance estimation. It should be taken into account that larger M leads to more accuracy at the expense of smaller susceptibility to dynamicity, and smaller M results in instability and inaccuracy [34]. The synergy of iterative and adaptive schemes in EKF improves the accuracy and robustness of the filter to nonlinearities and uncertainties of the system. The iterative approach recalculates the a posteriori estimate for a specified number of iterations. Afterwards, the adaptive method tunes the noise covariance matrices of the next time step based on previous measurements. Therefore, the filtering framework starts with an initial guess for F x , F y and F z . Then Equations (5) and (6) are integrated to provide k 1 , k 2 , τ. Curvatures and torsion are substituted in Equations (7) and (8). Obtaining the poses leads to specification of curvatures and torsions, TpX k q and ω se pX k q. Also, Y k is provided through the pose-measuring systems. The provided data are fed into the adopted estimating algorithm and therefore contact forces are estimated at each time step.

Simulation Results
In order to simulate the proposed model for tri-axial contact force sensing, an ablation catheter (Biosense Webster, Los Angeles, CA, USA) with the geometric and mechanical specifications as indicated in Table 2, is adopted. Simulations are conducted to illustrate and compare the accuracy and robustness of the developed contact force estimators (CFE) for in vivo catheterization, and then to represent the simultaneous contact force and stiffness estimations (CFSE) in required circumstances. At the beginning of the ablation procedure, the stiffness estimation is required while later during the intervention, for more computational efficiency, CFE is preferred.

Accuracy Analysis of CFE
In the first scenario, the accuracy of the CFE was tested. Tuning nonlinear Kalman filters is not a trivial task. To this end, measurement noise covariance matrix, R = diag(1 mm 2 , 1 mm 2 , 1 mm 2 , 0.0012 rad 2 , 0.0012 rad 2 , 0.0012 rad 2 , 1 mm´2, 1 mm´2, 1 mm´2) was obtained based on the accuracy of measuring device (Aurora ® EM Tracking System, NDI, Waterloo, ON, Canada). In the proposed model, Q can be tuned predicting the dynamics of states (which could vary significantly in practice). We chose Q " diagp9ˆ10´2 N 2 , 9ˆ10´2 N 2 , 9ˆ10´2 N 2 q. The human heart ordinarily beats 60-100 times per minute (1-1.7 Hz). Therefore, the designated observers should give reasonable response with sufficient sampling for such a dynamic environment. Figure 3 shows the tri-axial force exerted on the catheter tip including sinusoidal force with the abrupt change to the constant amount to entail regularity and irregularity in cardiac force exertion. Note that the steps represent sampling numbers in the corresponding figures. Figure 4 represents the corresponding errors of force estimation obtained by the EKF-based approaches when the filters are tuned initially. Table 3 gives the detailed statistical values of associated errors of the implemented methods. The results show that EKF does not converge to the ground truth value due to the high nonlinearity of the system as described in Section 3. However, iterative EKF (IEKF) and iterative adaptive EKF (IAEKF) show accurate performance, whereas IAEKF outperforms with the cost of higher computational time.   It should be noted that, when the distal shaft is straight, the axial contact force is ill-sensed. Figure 5 shows the insensible force along 3 d direction, when the distal shaft is straight. The true applied axial forces along with the estimated ones are represented in Figure 5, when the catheter moves with different initial curvature values, 1,0 , k varying in the range of 0-60 m −1 . However, since the distal shaft bends passing through the pulmonary vein, this limitation cannot be too restrictive.     It should be noted that, when the distal shaft is straight, the axial contact force is ill-sensed. Figure 5 shows the insensible force along 3 d direction, when the distal shaft is straight. The true applied axial forces along with the estimated ones are represented in Figure 5, when the catheter moves with different initial curvature values, 1,0 , k varying in the range of 0-60 m −1 . However, since the distal shaft bends passing through the pulmonary vein, this limitation cannot be too restrictive.   It should be noted that, when the distal shaft is straight, the axial contact force is ill-sensed. Figure 5 shows the insensible force along d 3 direction, when the distal shaft is straight. The true applied axial forces along with the estimated ones are represented in Figure 5, when the catheter moves with different initial curvature values, k 1,0 , varying in the range of 0-60 m´1. However, since the distal shaft bends passing through the pulmonary vein, this limitation cannot be too restrictive.

Convergence Analysis of CFE
The performance of a force estimator is usually evaluated according to the errors corresponding to the force estimation. Alternatively, convergence index (CI) can be calculated using the measurement, ˆr S r T k k k k   for EKF-based methods [31]. The smaller the value, the more stable the estimator is. Figure 6 illustrates the convergence of EKF-based approaches for the simulations performed in Section 4.1. It can be seen that IEKF and IAEKF have high convergence, where IAEKF excels in a smaller convergence index.

Convergence Analysis of CFE
The performance of a force estimator is usually evaluated according to the errors corresponding to the force estimation. Alternatively, convergence index (CI) can be calculated using the measurement, ρ k "r T k S krk for EKF-based methods [31]. The smaller the value, the more stable the estimator is. Figure 6 illustrates the convergence of EKF-based approaches for the simulations performed in Section 4.1. It can be seen that IEKF and IAEKF have high convergence, where IAEKF excels in a smaller convergence index.

Convergence Analysis of CFE
The performance of a force estimator is usually evaluated according to the errors corresponding to the force estimation. Alternatively, convergence index (CI) can be calculated using the measurement, ˆr S r T k k k k   for EKF-based methods [31]. The smaller the value, the more stable the estimator is. Figure 6 illustrates the convergence of EKF-based approaches for the simulations performed in Section 4.1. It can be seen that IEKF and IAEKF have high convergence, where IAEKF excels in a smaller convergence index.  show more robustness to changes of noise covariance matrix.

Sensitivity Analysis of CFE
(1) Measurement noise covariance matrix, R k : The convergence of KF-based approaches depends on the measurement noise covariance matrix. Five different measurement noise covariance matrices with diagonal elements R s P t0.02R k , 0.1R k , R k , 10R k , 50R k u with R k defined in Section 4.1, were adopted to study the sensitivity of the method to measurement noises. Table 4 compares the root mean square error (RMSE) for different methods with the applied force along x direction, F x " 0.1sinp2πtq, with 15 samples per second. As it can be seen, IAEKF and IEKF show more robustness to changes of noise covariance matrix. (2) Initial state estimate,X 0 : The cardiac ablation force varies within the range of (0.1-0.4 N) [1].
Therefore, our initial estimates cannot be too far from this range for cardiac-intervention applications. Two sets of initial values, one with 20% (0.2 F max ) and the other one with 50% (0.5 F max ) of the maximum true values, were used as initial estimates. Table 5 compares the corresponding errors of two different initializations for the applied force of F x " 0.1sinp2πtq.
The sensitivity analyses demonstrate that IEKF and IAEKF have the highest robustness to the filtering uncertainties that might occur during the ablation.

Contact Force and Stiffness Estimation
In previous studies, stiffness was assumed to be as an invariant parameter of the system [15,17]. However, this parameter might change as a function of structural specifications of the catheter and different environmental conditions, e.g., environmental temperature [22]. For this purpose, CFSE was tested. The process covariance matrix was assumed to be Q " diagp9ˆ10´4 N 2 , 9ˆ10´4 N 2 , 91 0´4 N 2 , 10´8 N 2 m 4 , 10´8 N 2 m 4 q, and R was defined equal to that of CFE in Section 4.1. With the tip contact force represented in Figure 3, the error of force and stiffness estimations are indicated in Figure 7a,b, respectively. The results demonstrate that IEKF and IAEKF provide estimates of contact force with the maximum average error of´0.0016˘0.0036 and 0.0008˘0.0054, respectively, and stiffness estimates with´1.59ˆ10´4 and 1.28ˆ10´5, respectively. The results for EKF-based stiffness estimation due to its divergence have been eliminated. The increased inaccuracy of force estimation is mostly due to high error at the first instants of estimation. Although there is an increased error of contact force estimation compared to that of CFE, the proper choice of measurement vector yielded to convergent estimations for both IEKF and IAEKF. However, depending on the parameters to be estimated, the number of measurement points could be increased.
In the end, the execution times of the proposed approaches are compared. Table 6 illustrates the mean value of computational time of the methods for the simulations. The simulations are performed in MATLAB 2013a ® with Core i3 2.2 GHz, 4GB RAM laptop. The results show that IEKF, IAEKF take longer execution time, mostly due to the iterative nature of the schemes. Numerical Jacobian calculation leads to higher computational time that can be compensated by exploiting faster computers, and running the codes in C/C++.

Experimental Validation
Experiments were designed to compare and verify the performance of the proposed methods. To obtain the stiffness of the catheter, tip forces were applied to the distal shaft and the shape of it was obtained through the measurements of equally spaced point positions 1 cm apart on the distal shaft using CMM machine (Mitutoyo AE112, resolution: 0.001 mm). In our analysis, we considered constant bending stiffness for distal shaft in both directions and The optimization was performed using Nelder-Mead simplex method within MATLAB. The results of identification for the ThermoCool SF ablation catheter are presented in Table 2.

Experiment I
In the first experiment, calibration weights were applied to the catheter tip. Figure 8a represents the experimental setup provided for tri-axial load exertion simulating the forces applied in ablation procedure. The setup consisted of a thread, connected to the tip of the distal shaft, passing through the frictionless pulley with the calibrating weights hung on the other end of the thread. In order to apply spatial forces with the desired direction, the pulley was mounted on the framework with two moving links, as depicted in Figure 8a.

Experimental Validation
Experiments were designed to compare and verify the performance of the proposed methods. To obtain the stiffness of the catheter, tip forces were applied to the distal shaft and the shape of it was obtained through the measurements of equally spaced point positions 1 cm apart on the distal shaft using CMM machine (Mitutoyo AE112, resolution: 0.001 mm). In our analysis, we considered constant bending stiffness for distal shaft in both directions and GJ s " 0.3846 EpI 1`I2 q with v p " 0.3. The identification procedure is accomplished by conducting n e " 10 experiments and solving nonlinear optimization problems for two unknown variables tEI 1 , EI 2 u by minimizing the distance between measured and model predicted data points as follows: min pJ cost q " min¨n The optimization was performed using Nelder-Mead simplex method within MATLAB. The results of identification for the ThermoCool SF ablation catheter are presented in Table 2.

Experiment I
In the first experiment, calibration weights were applied to the catheter tip. Figure 8a represents the experimental setup provided for tri-axial load exertion simulating the forces applied in ablation procedure. The setup consisted of a thread, connected to the tip of the distal shaft, passing through the frictionless pulley with the calibrating weights hung on the other end of the thread. In order to apply spatial forces with the desired direction, the pulley was mounted on the framework with two moving links, as depicted in Figure 8a.
The system in Figure 8b was used to examine the applicability of force estimation approaches. Aurora electromagnetic (EM) tracking system (NDI, Waterloo, ON) consisting of EM field generator and EM tracking sensors provided the position and orientation (pose) of its sensors with the accuracy of 0.48 mm and 0.3˝. The curvatures and torsion were obtained through calculating the changes of orientation over arc length in two adjacent end segments of the distal shaft with 4.7 mm distance from one another. Sensor coils were incorporated for measuring the pose of the distal tip and the marked points on the distal shaft as well. Three-dimensional forces were estimated using the measurements and the proposed KF-based approaches.
Sensors 2016, 16,990 14 of 20 The system in Figure 8b was used to examine the applicability of force estimation approaches. Aurora electromagnetic (EM) tracking system (NDI, Waterloo, ON) consisting of EM field generator and EM tracking sensors provided the position and orientation (pose) of its sensors with the accuracy of 0.48 mm and 0.3°. The curvatures and torsion were obtained through calculating the changes of orientation over arc length in two adjacent end segments of the distal shaft with 4.7 mm distance from one another. Sensor coils were incorporated for measuring the pose of the distal tip and the marked points on the distal shaft as well. Three-dimensional forces were estimated using the measurements and the proposed KF-based approaches. The catheter distal shaft was bent into two different initial curvatures. The curvature was calculated through marked equispaced point position measurements with 1 cm distance and fitting the optimal constant curvature arc by means of Nelder-Mead simplex method using MATLAB's fminsearch function [35]. In the first scenario of experiment I, planar forces were applied in different directions with the weight 13.03 grams (resembling moderate ablation force of 0.13 N) hung to the pulley and rotating the second link. In the second scenario of experiment I, the spatial forces were applied with the same weight as in previous experiment, but rotating the two links resembling the sinusoidal functions in 3D. The applied forces along the local axes of the distal shaft tip were obtained via pose measurements and obtaining the projected weight forces along the local axes. Figure 9 shows the true values of planar and spatial forces applied to the end of the distal shaft with initial curvatures of respectively. Each experiment was repeated three times for verifying the repeatability of the introduced methods. Figure 10a,b represent the corresponding errors of estimation for aforementioned contact forces performed by EKF, IEKF, and IAEKF. For better comparison of the results, the average error values and corresponding standard deviation of estimations are summarized in Table 7. The entries of the measurement covariance matrix were taken to be 1000 times that of the determined R k to account for the measurement uncertainty due to high sensitivity of the EM tracking system to trivial distortions. The results show that IAEKF has the highest accuracy with the maximum average error value along axial direction with the precision of 0.001 ± 0.0041 N. However, the force estimations provided by IEKF and IAEKF are comparable with the existing force sensors [11] and estimators [18]. The catheter distal shaft was bent into two different initial curvatures. The curvature was calculated through marked equispaced point position measurements with 1 cm distance and fitting the optimal constant curvature arc by means of Nelder-Mead simplex method using MATLAB's fminsearch function [35]. In the first scenario of experiment I, planar forces were applied in different directions with the weight 13.03 grams (resembling moderate ablation force of 0.13 N) hung to the pulley and rotating the second link. In the second scenario of experiment I, the spatial forces were applied with the same weight as in previous experiment, but rotating the two links resembling the sinusoidal functions in 3D. The applied forces along the local axes of the distal shaft tip were obtained via pose measurements and obtaining the projected weight forces along the local axes. Figure 9 shows the true values of planar and spatial forces applied to the end of the distal shaft with initial curvatures of k 1,0 " 39.99 m´1 and k 1,0 " 24.57 m´1, respectively. Each experiment was repeated three times for verifying the repeatability of the introduced methods. Figure 10a,b represent the corresponding errors of estimation for aforementioned contact forces performed by EKF, IEKF, and IAEKF. For better comparison of the results, the average error values and corresponding standard deviation of estimations are summarized in Table 7. The entries of the measurement covariance matrix were taken to be 1000 times that of the determined R k to account for the measurement uncertainty due to high sensitivity of the EM tracking system to trivial distortions. The results show that IAEKF has the highest accuracy with the maximum average error value along axial direction with the precision of 0.001˘0.0041 N. However, the force estimations provided by IEKF and IAEKF are comparable with the existing force sensors [11] and estimators [18].

Experiment II
In the second set of experiments, CFSE and CFE were evaluated in a simulated environment of intra-cardiac catheterization. For CFE, the setup in Figure 11a was prepared. The water circulation pump (EHIEM, Deizisau, Germany) circulated the saline water in the bucket mimicking the blood flow in intra-cardiac chambers. The force sensor (Nano 17, ATI, Apex, NC, USA) was attached to the bucket wall to interact with the catheter tip while measuring the tri-axial contact forces. The catheter was bent with initial curvature of

Experiment II
In the second set of experiments, CFSE and CFE were evaluated in a simulated environment of intra-cardiac catheterization. For CFE, the setup in Figure 11a was prepared. The water circulation pump (EHIEM, Deizisau, Germany) circulated the saline water in the bucket mimicking the blood flow in intra-cardiac chambers. The force sensor (Nano 17, ATI, Apex, NC, USA) was attached to the bucket wall to interact with the catheter tip while measuring the tri-axial contact forces. The catheter was bent with initial curvature of

Experiment II
In the second set of experiments, CFSE and CFE were evaluated in a simulated environment of intra-cardiac catheterization. For CFE, the setup in Figure 11a was prepared. The water circulation pump (EHIEM, Deizisau, Germany) circulated the saline water in the bucket mimicking the blood flow in intra-cardiac chambers. The force sensor (Nano 17, ATI, Apex, NC, USA) was attached to the bucket wall to interact with the catheter tip while measuring the tri-axial contact forces. The catheter was bent with initial curvature of k 1,0 " 19.41 m´1. Two Aurora sensors were attached to the catheter distal shaft to provide tip pose, curvatures and torsion of the second-last segment required for force sensing as illustrated in Figure 11b. The sensors were placed 4.7 mm apart and calibrated to measure the desired poses of the two end segments.
Sensors 2016, 16,990 16 of 20 sensing as illustrated in Figure 11b. The sensors were placed 4.7 mm apart and calibrated to measure the desired poses of the two end segments. In this experiment, the bucket with the force sensor was moved forward manually toward the distal shaft base, and the force estimation was performed for the sampled data points. Figure 12 illustrates the applied and estimated spatial contact forces expressed in the local distal shaft frame coordinates. In the second set of experiments, the measurement covariance matrix was assumed to be equal to that of Experiment I, except the last three diagonal elements when it was assumed to be twice those of the covariance matrix due to uncertainties imposed by circulating water over the curvature and torsion measurements of the second-last segment. In addition, the estimations for the sampled data points were repeated five times. The results show the consistency with the previous experiments and simulations, where the IAEKF estimator has the maximum average error of 0.0028 ± 0.0019 N along the axial direction. The average errors of estimations are summarized in Table 7. In the second scenario, the setup in Figure 13 was prepared for CFSE. To simulate the beating heart as the dynamic interacting contact load of the distal shaft, two DC servomotors (The Lego Group, Billund, Denmark), a bucket accommodated on top of the motors, and the force sensor attached to the bucket wall, as represented in Figure 13a,b, were incorporated. The distal shaft tip was in contact with the force sensor and the motors were connected together and were moving the In this experiment, the bucket with the force sensor was moved forward manually toward the distal shaft base, and the force estimation was performed for the sampled data points. Figure 12 illustrates the applied and estimated spatial contact forces expressed in the local distal shaft frame coordinates. In the second set of experiments, the measurement covariance matrix was assumed to be equal to that of Experiment I, except the last three diagonal elements when it was assumed to be twice those of the covariance matrix due to uncertainties imposed by circulating water over the curvature and torsion measurements of the second-last segment. In addition, the estimations for the sampled data points were repeated five times. The results show the consistency with the previous experiments and simulations, where the IAEKF estimator has the maximum average error of 0.0028˘0.0019 N along the axial direction. The average errors of estimations are summarized in Table 7.
Sensors 2016, 16, 990 16 of 20 sensing as illustrated in Figure 11b. The sensors were placed 4.7 mm apart and calibrated to measure the desired poses of the two end segments. In this experiment, the bucket with the force sensor was moved forward manually toward the distal shaft base, and the force estimation was performed for the sampled data points. Figure 12 illustrates the applied and estimated spatial contact forces expressed in the local distal shaft frame coordinates. In the second set of experiments, the measurement covariance matrix was assumed to be equal to that of Experiment I, except the last three diagonal elements when it was assumed to be twice those of the covariance matrix due to uncertainties imposed by circulating water over the curvature and torsion measurements of the second-last segment. In addition, the estimations for the sampled data points were repeated five times. The results show the consistency with the previous experiments and simulations, where the IAEKF estimator has the maximum average error of 0.0028 ± 0.0019 N along the axial direction. The average errors of estimations are summarized in Table 7. In the second scenario, the setup in Figure 13 was prepared for CFSE. To simulate the beating heart as the dynamic interacting contact load of the distal shaft, two DC servomotors (The Lego Group, Billund, Denmark), a bucket accommodated on top of the motors, and the force sensor attached to the bucket wall, as represented in Figure 13a,b, were incorporated. The distal shaft tip was in contact with the force sensor and the motors were connected together and were moving the In the second scenario, the setup in Figure 13 was prepared for CFSE. To simulate the beating heart as the dynamic interacting contact load of the distal shaft, two DC servomotors (The Lego Group, Billund, Denmark), a bucket accommodated on top of the motors, and the force sensor attached to the bucket wall, as represented in Figure 13a,b, were incorporated. The distal shaft tip was in contact with the force sensor and the motors were connected together and were moving the bucket. The bucket was translated with a frequency and sampling time of approximately 0.21 Hz and 0.3 s, respectively, to comply with the computational requirements of CFSE framework. However, the computational costs easily can be reduced as mentioned in Section 4. In cardiac catheterization, the blood flow has the most significant effect in making disturbances in both sensor readings and process model. In this scenario, water was circulated with the other DC servomotor as illustrated in Figure 13a. Force and EM sensors were calibrated and attached to the bucket wall and distal shaft, respectively, as in the first scenario of Experiment II. bucket. The bucket was translated with a frequency and sampling time of approximately 0.21 Hz and 0.3 s, respectively, to comply with the computational requirements of CFSE framework. However, the computational costs easily can be reduced as mentioned in Section 4. In cardiac catheterization, the blood flow has the most significant effect in making disturbances in both sensor readings and process model. In this scenario, water was circulated with the other DC servomotor as illustrated in Figure 13a. Force and EM sensors were calibrated and attached to the bucket wall and distal shaft, respectively, as in the first scenario of Experiment II.
(a) (b) Figure 13. Experimental apparatus mimicking an intra-cardiac chamber for CFSE: (a) Close-up view of the setup consisting of EM tracking system, a motor for fluid circulation, ablation catheter and the force sensor; (b) Two dc servomotors translating the aforementioned assembly for making the dynamic interacting load. Figure 14a,b represent the results of contact force and stiffness estimation of distal shaft interacting with the moving force sensor. There is no direct way for evaluation of the accuracy of stiffness estimation in different operational conditions. However, it can be assessed indirectly from the accuracy of force estimation. The maximum average error of force estimation for IEKF and IAEKF are 0.0138 ± 0.0050 N and 0.0035 ± 0.0012 N, respectively, where the IAEKF estimator due to its adaptability has the highest accuracy and robustness to the uncertainties. Despite shortcomings of the experiment (i.e., ignoring cardiac tissue elasticity, catheter's pathway from pulmonary vein to the atrial and compliance of blood vessels connected to the heart), it provides good basis for proof of concept. It should also be noted that this is the first study considering simultaneous force and stiffness estimation.    Figure 14a,b represent the results of contact force and stiffness estimation of distal shaft interacting with the moving force sensor. There is no direct way for evaluation of the accuracy of stiffness estimation in different operational conditions. However, it can be assessed indirectly from the accuracy of force estimation. The maximum average error of force estimation for IEKF and IAEKF are 0.0138˘0.0050 N and 0.0035˘0.0012 N, respectively, where the IAEKF estimator due to its adaptability has the highest accuracy and robustness to the uncertainties. Despite shortcomings of the experiment (i.e., ignoring cardiac tissue elasticity, catheter's pathway from pulmonary vein to the atrial and compliance of blood vessels connected to the heart), it provides good basis for proof of concept. It should also be noted that this is the first study considering simultaneous force and stiffness estimation. bucket. The bucket was translated with a frequency and sampling time of approximately 0.21 Hz and 0.3 s, respectively, to comply with the computational requirements of CFSE framework. However, the computational costs easily can be reduced as mentioned in Section 4. In cardiac catheterization, the blood flow has the most significant effect in making disturbances in both sensor readings and process model. In this scenario, water was circulated with the other DC servomotor as illustrated in Figure 13a. Force and EM sensors were calibrated and attached to the bucket wall and distal shaft, respectively, as in the first scenario of Experiment II.
(a) (b) Figure 13. Experimental apparatus mimicking an intra-cardiac chamber for CFSE: (a) Close-up view of the setup consisting of EM tracking system, a motor for fluid circulation, ablation catheter and the force sensor; (b) Two dc servomotors translating the aforementioned assembly for making the dynamic interacting load. Figure 14a,b represent the results of contact force and stiffness estimation of distal shaft interacting with the moving force sensor. There is no direct way for evaluation of the accuracy of stiffness estimation in different operational conditions. However, it can be assessed indirectly from the accuracy of force estimation. The maximum average error of force estimation for IEKF and IAEKF are 0.0138 ± 0.0050 N and 0.0035 ± 0.0012 N, respectively, where the IAEKF estimator due to its adaptability has the highest accuracy and robustness to the uncertainties. Despite shortcomings of the experiment (i.e., ignoring cardiac tissue elasticity, catheter's pathway from pulmonary vein to the atrial and compliance of blood vessels connected to the heart), it provides good basis for proof of concept. It should also be noted that this is the first study considering simultaneous force and stiffness estimation.  The primary sources of errors in the experiments are the pose measuring system and consequently deformation calculation inaccuracy, and discrepancy of evaluated initial curvatures from the true values. The former one can be tackled with the exploitation of more accurate measuring system and imaging modalities, while it can relax the inaccuracy of distal shaft initial curvature estimation as well. In addition, incorporating redundant pose measurements along the distal shaft can lead to more precise initial curvature attribution. In clinical applications, image-based pose measurements would be feasible. Modeling uncertainty is the other source of error. The model can be elaborated by considering the contact forces and tendon loadings acting along the distal shaft. However, the computational time and accuracy tradeoff should be taken into account.

Discussion
In intra-cardiac catheterization procedures' force measurements are very helpful for the success of both manual and robot-assisted catheterization. A new framework for tri-axial shape-based contact force and stiffness estimation of tendon-driven catheters was introduced in this paper. The introduced estimators show the potential of Kalman filter (KF)-and observer-based methods as the alternative for integrating often costly and inconvenient catheter tip force sensors. In particular, readily available images during intervention can be used to estimate the curvatures and position of required points on the catheter for the proposed method. It has been demonstrated that contact forces can be observed using catheter deformation [18].
In this work, nonlinear Kalman filters utilizing tip pose and curvature measurements were formulated for force and stiffness estimation. The shape-based force estimators were developed based on the accurate Krichhoff rod model of the tendon-driven catheter. To cope with model and measurement nonlinearities and uncertainties, iterative and adaptive structures were integrated into EKF formulation to provide reliable force and stiffness estimators. The IAEKF approach increases the accuracy and robustness of the estimation considerably at the expense of more computational time. Implementation of iterative approaches requires tuning of iteration numbers, since over-iterating results in higher computational time. The window size of adaptive filters should be adjusted as well, since its improper size would result in the instability of the system. The simulation and experimental results show the potential application of the proposed method in clinical experiments. For future work, we intend to perform in vivo tests. For achieving more accurate and redundant pose measurements, we plan to use medical imaging modalities along with EM tracking sensor measurements. The filtering framework is expected to be elaborated upon in the presence of point and distributed forces along the distal shaft taking into account changing contact forces as a result of varying heart beat rate in cardiac arrhythmia treatment. Furthermore, we plan to incorporate measurements of tissue deformation in the case of zero curvature of the distal tip.