Front. Robot. AI, 30 July 2021
Sec. Biomedical Robotics

FBG-Based Estimation of External Forces Along Flexible Instrument Bodies

  • 1Robot Assisted Surgery (RAS), Department of Mechanical Engineering, KU Leuven University, Leuven, Belgium
  • 2FBGS International NV, Geel, Belgium

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 De Greef et al. (2009); Burgner-Kahrs et al. (2015); Heunis et al. (2018). 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 Balami et al. (2018); Lv et al. (2011); Kavic and Basson (2001). As a consequence, the demand for sensorized instruments relaying crucial real-time information such as: position, shape, and interaction force is ever increasing Shi et al. (2017); Wu et al. (2021); Galloway et al. (2019); Trejos et al. (2010); Haidegger et al. (2009).

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 (Jung et al., 2006; Colome et al., 2013; Daly and Wang, 2014; Hu and Xiong, 2018), and more recently, on soft/continuum robots (Yasin and Simaan, 2020; Sadati et al., 2020; Xu and Simaan, 2010; Gao et al., 2020). 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 (Rucker and Webster, 2011). 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 (Back et al., 2015). Hooshiar et al. estimate tip forces by employing Bezier spline shape approximations and an inverse Cosserat rod model (Hooshiar et al., 2020). Sadati et al. employed force sensor readings at the base of a continuum appendage and the Cosserat rod model for tip load estimation (Sadati et al., 2020). Further examples of previous works on tip force estimation using kinematic models can be found in (Khoshnam et al., 2015; Khan et al., 2017).

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 (Qiao et al., 2019). 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 (Qiao et al., 2021). 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 (Aloi and Rucker, 2019). Heunis et al. (Heunis et al., 2019) and Venkiteswaran et al. (Venkiteswaran et al., 2019) 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 (Ryu et al., 2020). 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 (Massari et al., 2020). 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 λB 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 ΔT. Considering small temperature shifts, the change in Bragg wavelength can be expressed as:


where λB0 is the grating’s unstrained Bragg wavelength, Δϵ is the change in strain with respect to the unstrained state, and Sϵ and ST 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 (ΔλB,1). For the outer cores, the temperature contribution STΔT is thus known from the central core, and Eq. 1 can be rewritten as:


where Δϵi is the change in strain in the ith core. The relationship between the strain in each core ϵi 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 θi is the angle of the ith. 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 θb, where θb 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 θb by defining an apparent curvature vector κappMoore and Rogge (2012); Al-Ahmad et al. (2020):

κ=2 |κapp|N,(5)

where i^ and j^ 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 θb in Eq. 6 is obtained using the atan2 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 κ[j] and bend angles θb[j] along the fiber’s length with j[1,n].


FIGURE 1. (A) cross-sectional view of a four core MCF where λi represents the wavelength in the ith core, r is the distance between the center of the ith core and the fiber’s central axis, and θb is the angle of the bending plane with respect to the x-axis which intersects the second core, (B) isometric view of a segment with four cores and two FBG sets, i.e. four FBGs at a given cross-section, separated by a center-to-center distance lz.

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. A deformed elastic rod parametrized by arc length s subjected to a set of point forces F1 and F2, and distributed forces f1 and f2.

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 Antman (2005):


where n and m are [3×1] internal force and internal moment vectors, respectively, u^ and v^ are the [3×3] angular and linear strain matrices in skew-symmetric form, respectively, and f and l are the [3×1] 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 v and u are the [3×1] linear and angular strain vectors, respectively, Kv1 is the [3×3] shear and axial stiffness matrix, Ku1 is the [3×3] bending and torsional stiffness matrix, and v0 and u0 are the [3×1] reference linear and angular strain vectors, respectively. Note that the angular strain u in Eq. 10 corresponds to:


where κ is the curvature given in Eq. 5 and θb 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. v=v0=[0,0,1]T, such that axial forces are disregarded

• there is no, or negligible torsion, i.e. u[3]=0

• 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 n and internal moment m at the rod’s “boundaries”, i.e. at the base of the rod s=0, and at the tip of the rod s=L, 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 s=0, it is known that the internal force n and internal moment m are both zero at s=L 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 s=L with boundary conditions n=0 and m=0 (or n=Fext and m=Mext), and solving backwards towards s=0.

Notice that the set of equilibrium equations do not contain terms for externally applied point forces F nor point moments M. 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 FL/2 and a point moment ML/3 at locations s=L/2 and s=L/3, respectively. The solution can then be found by dividing the problem into three integration sets, starting with n=0 and m=0 at s=L, then n=n++FL/2 at s=L/2, and finally m=m++ML/3 at s=L/3. Here, the - and + indicates the value of n or m just before the given location (towards s=0) and just after the given location (towards s=L), 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 m can be computed in Eq. 10 given knowledge of the angular strain vector u. Notice that the stiffness matrix Ku and reference angular strain u0 are constants, and that u 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 n (given that l=0 was assumed). Generally, the resulting force profile of n 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 Fext could be found at the break points (see Figure 3), i.e. points of discontinuity, in the internal force profile such that Fext=nn+. 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. Example plot illustrating the differences between the theoretical, measured, and estimated curvature profiles.

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 n=1 and would result in a linear curvature profile, n=1. Similarly, a distributed force with constant magnitude has a degree n=0 and would result in a quadratic curvature profile n=2. Hence, to overcome the noisy curvature measurements problem stated previously, the curvature profile can be segmented into a series of 1st order polynomials for point loads, or a series of 1st and 2nd 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 Lovrić et al. (2014) 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 Qiao et al. (2019). 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 Lovrić et al. (2014).

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 [1...Nsegs,max], 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 ϵk,thresh.


TABLE 1. 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 Nsegs. 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 Fmin (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 Killick et al. (2012). This time however, instead of defining the number of segments to find the break points, a minimum slope change Δdκ/dsmin 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 segsOut are used in the constitutive Eq. 10 to compute the internal moment m. 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 Nact. If it is strictly larger, it will become clear that the magnitudes of those forces additional to Nact will be equal to zero, or have a very small value. Let us now consider a function g(χ) that requires an input χ and returns an output κ such that:


where Fxi and Fyi are the magnitudes of the x and y components of the ith external point force relative to the rod’s body frame, respectively. Li is the location of the external point force expressed as arc length s=Li. κxj and κyj are the magnitudes of the x and y components of the jth curvature measurement along the arc length, and M the total number of curvature measurements over the total length. The function g(χ) 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 g(χ).

3.3.1 Least-Squares Optimization

An optimization algorithm typically finds a set of optimization parameters based on minimizing a desired cost function Θ(χ,κmeas). In this case, the cost function is defined as:


where κmeas is the measured curvature vector using the MCF. Hence, the goal is to minimize the error between the measured curvature κmeas and the estimated curvature g(χ)=κ, 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 Le et al. (2017).

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 Wan and Van Der Merwe (2000). 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 xk are the states at time step k such that xk=χ, h(xk) is the measurement function such that h(xk)=g(xk), zk is the measurement at time step k, and wkN(0,Qk) and vkN(0,Rk) are the process and measurement noise, respectively. Qk and Rk are the process and measurment noise covariance matrices at time step k, and are found by manual tuning. Note that the states xk, 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 Wan and Van Der Merwe (2000), 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 lz 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. Force estimation experimental setup comprising: ① electromagnetic (EM) field generator, ② electromagnetic tracking (EM) sensors, ③ ATI Nano17 six DoF force sensor, ④ tensioned cable to transmit a force to the Nitinol rod, ⑤ weights suspended by gravity, ⑥ Nitinol rod with an MCF embedded internally, ⑦ optical fanout and interrogator systems, ⑧ monitor to display wavelength spectra. Forces F1 and F3 are applied by suspending known weights, and force F2 is applied along an adjustable direction by pulling on the ATI Nano17 six DoF force sensor which is attached to a cable that encircles the Nitinol rod.

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 Ftot,ATI=Fx,ATI2+Fy,ATI2+Fz,ATI2=F2, where F2 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 F2 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 s=0, and the second is attached at the tip of the rod at s=L. 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 Al-Ahmad et al. (2020) for the method used to reconstruct the rod’s shape]. The rigid transformation matrix TrodEM 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 F2, 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: ϵκ,thresh=0.05, Nsegs,max=5, Fmin=0.2 N, Nsegs=5, and Δdκ/dsmin=0.45 m−1. Note that ϵκ,thresh is a normalized value. Furthermore, N=4, Rk = diag [repmat (0.01 m−1, 1, 2M)], and Qk = diag [repmat (0.1 N, 0.1 N, 0.1 m, 1, N)] where the function repmat (matIn, Nrow, Ncol) returns a matrix that is a repetition of the input matrix matIn for Nrow rows and Ncol columns, and the function diag (matIn) returns a diagonal matrix based on the elements of matIn.

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 Fmin (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. Summary of force estimation results illustrating estimated force numbers and errors on magnitudes and locations. E(||||) 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 θb 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 Fmin, 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. Example result for the case where three external point forces are applied onto the rod. (A) profile of the x-component of the curvature, (B) profile of the y-component of the curvature, (C) deformed rod showing a comparison between the ground truth and estimated shape and applied forces.

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 fxi and fyi are the magnitudes of the x and y components of the ith external distributed force, respectively, and Lsi and Lei are the start and end locations of the distributed force along the arc length. The output curvature vector κ and the forward kinematics function g(χ) remain the same. The only difference is that distributed forces f 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. UKF simulation with two distributed forces applied onto a rod. The magnitudes and locations are varying sinusoidally over time (gt = ground truth, est = estimated).


TABLE 3. 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.

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.


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.


Al-Ahmad, O., Ourak, M., Van Roosbroeck, J., Vlekken, J., and Vander Poorten, E. B. (2020). Improved FBG-Based Shape Sensing Methods for Vascular Catheterization Treatment. IEEE Robot. Autom. Lett. 5, 1. doi:10.1109/LRA.2020.3003291

CrossRef Full Text | Google Scholar

Aloi, V. A., and Rucker, D. C. (2019). “Estimating Loads along Elastic Rods,” in IEEE International Conference on Robotics and Automation, Montreal, QC, Canada, May 20–24, 2019, 2867–2873. doi:10.1109/ICRA.2019.8794301

CrossRef Full Text | Google Scholar

Antman, S. S. (2005). Nonlinear Problems of Elasticity. Appl. Math. Sci. 107, 1–835. doi:10.1007/978-1-4757-4147-6_1

CrossRef Full Text | Google Scholar

Back, J., Manwell, T., Karim, R., Rhode, K., Althoefer, K., and Liu, H. (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. doi:10.1109/IROS.2015.7353647

CrossRef Full Text | Google Scholar

Balami, J. S., White, P. M., McMeekin, P. J., Ford, G. A., and Buchan, A. M. (2018). Complications of Endovascular Treatment for Acute Ischemic Stroke: Prevention and Management. Int. J. Stroke 13, 348–361. doi:10.1177/1747493017743051

PubMed Abstract | CrossRef Full Text | Google Scholar

Burgner-Kahrs, J., Rucker, D. C., and Choset, H. (2015). Continuum Robots for Medical Applications: A Survey. IEEE Trans. Robot. 31, 1261–1280. doi:10.1109/TRO.2015.2489500

CrossRef Full Text | Google Scholar

Colome, A., Pardo, D., Alenya, G., and Torras, C. (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. doi:10.1109/ICRA.2013.6631072

CrossRef Full Text | Google Scholar

Daly, J. M., and Wang, D. 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. doi:10.1109/TCST.2013.2242329

CrossRef Full Text | Google Scholar

De Greef, A., Lambert, P., and Delchambre, A. (2009). Towards Flexible Medical Instruments: Review of Flexible Fluidic Actuators. Precision Eng. 33, 311–321. doi:10.1016/j.precisioneng.2008.10.004

CrossRef Full Text | Google Scholar

Galloway, K. C., Chen, Y., Templeton, E., Rife, B., Godage, I. S., and Barth, E. J. (2019). Fiber Optic Shape Sensing for Soft Robotics. Soft Robot. 6, 671–684. doi:10.1089/soro.2018.0131

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, A., Liu, N., Shen, M., E.M.K. Abdelaziz, M., Temelkuran, B., and Yang, G.-Z. (2020). Laser-Profiled Continuum Robot with Integrated Tension Sensing for Simultaneous Shape and Tip Force Estimation. Soft Robot. 7, 421–443. doi:10.1089/soro.2019.0051

PubMed Abstract | CrossRef Full Text | Google Scholar

Haidegger, T., Benyó, B., Kovács, L., and Benyó, Z. (2009). Force Sensing and Force Control for Surgical Robots. IFAC Proc. Volumes 42, 401–406. doi:10.3182/20090812-3-DK-2006.0035

CrossRef Full Text | Google Scholar

Hasanzadeh, S., and Janabi-Sharifi, F. (2015). Model-based Force Estimation for Intra-cardiac Catheters. IEEE ASME Trans. Mechatron. 21, 1. doi:10.1109/TMECH.2015.2453122

CrossRef Full Text | Google Scholar

Heunis, C. M., Belfiore, V., Vendittelli, M., and Misra, S. (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. doi:10.1109/IROS40897.2019.8968448

CrossRef Full Text | Google Scholar

Heunis, C., Sikorski, J., and Misra, S. (2018). Flexible Instruments for Endovascular Interventions: Improved Magnetic Steering, Actuation, and Image-Guided Surgical Instruments. IEEE Robot. Automat. Mag. 25, 71–82. doi:10.1109/MRA.2017.2787784

CrossRef Full Text | Google Scholar

Hooshiar, A., Sayadi, A., Jolaei, M., and Dargahi, J. (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. doi:10.1109/BIA50171.2020.9244512

CrossRef Full Text | Google Scholar

Hu, J., and Xiong, R. (2018). Contact Force Estimation for Robot Manipulator Using Semiparametric Model and Disturbance Kalman Filter. IEEE Trans. Ind. Electron. 65, 3365–3375. doi:10.1109/TIE.2017.2748056

CrossRef Full Text | Google Scholar

Jung, J., Lee, J., and Huh, K. (2006). Robust Contact Force Estimation for Robot Manipulators in Three-Dimensional Space. Proc. Inst. Mech. Eng. C. J. Mech. Eng. Sci. 220, 1317–1327. doi:10.1243/09544062C09005

CrossRef Full Text | Google Scholar

Kavic, S. M., and Basson, M. D. (2001). Complications of Endoscopy. Am. J. Surg. 181, 319–332. doi:10.1016/S0002-9610(01)00589-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, F., Roesthuis, R. J., and Misra, S. (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. doi:10.1109/IROS.2017.8206073

CrossRef Full Text | Google Scholar

Khoshnam, M., Skanes, A. C., and Patel, R. V. (2015). Modeling and Estimation of Tip Contact Force for Steerable Ablation Catheters. IEEE Trans. Biomed. Eng. 62, 1404–1415. doi:10.1109/TBME.2015.2389615

PubMed Abstract | CrossRef Full Text | Google Scholar

Killick, R., Fearnhead, P., and Eckley, I. A. (2012). Optimal Detection of Changepoints with a Linear Computational Cost. J. Am. Stat. Assoc. 107, 1590–1598. doi:10.1080/01621459.2012.737745

CrossRef Full Text | Google Scholar

Le, T. M., Fatahi, B., Khabbaz, H., and Sun, W. (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. doi:10.1016/j.apm.2016.08.034

CrossRef Full Text | Google Scholar

Lovrić, M., Milanović, M., and Stamenković, M. (2014). Algoritmic Methods for Segmentation of Time Series: An Overview. J. Contemp. Econ. Bus. Issues 1, 31–53.

Google Scholar

Lv, X., Wu, Z., Jiang, C., Li, Y., Yang, X., and Zhang, Y. (2011). Complication Risk of Endovascular Embolization for Cerebral Arteriovenous Malformation. Eur. J. Radiol. 80, 776–779. doi:10.1016/j.ejrad.2010.09.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Massari, L., Schena, E., Massaroni, C., Saccomandi, P., Menciassi, A., Sinibaldi, E., 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. doi:10.1089/soro.2018.0172

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, J. P., and Rogge, M. D. (2012). Shape Sensing Using Multi-Core Fiber Optic cable and Parametric Curve Solutions. Opt. Express 20, 2967–2973. doi:10.1364/oe.20.002967

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiao, Q., Borghesan, G., De Schutter, J., and Vander Poorten, E. (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.

Google Scholar

Qiao, Q., Willems, D., Borghesan, G., Ourak, M., De Schutter, J., and Poorten, E. 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. doi:10.1109/ICAR46387.2019.8981592

CrossRef Full Text | Google Scholar

Rucker, D. C., and Webster, R. 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. doi:10.1109/IROS.2011.6048202

CrossRef Full Text | Google Scholar

Ryu, H. T., Woo, J., So, B. R., and Yi, B. J. (2020). Shape and Contact Force Estimation of Inserted Flexible Medical Device. Int. J. Control Autom. Syst. 18, 163–174. doi:10.1007/s12555-019-0237-8

CrossRef Full Text | Google Scholar

Sadati, S. M., Shiva, A., Herzig, N., Rucker, C. D., Hauser, H., Walker, I. 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. doi:10.1109/LRA.2020.2972790

CrossRef Full Text | Google Scholar

Shi, C., Luo, X., Qi, P., Li, T., Song, S., Najdovski, Z., et al. (2017). Shape Sensing Techniques for Continuum Robots in Minimally Invasive Surgery: A Survey. IEEE Trans. Biomed. Eng. 64, 1665–1678. doi:10.1109/TBME.2016.2622361

PubMed Abstract | CrossRef Full Text | Google Scholar

Trejos, A. L., Patel, R. V., and Naish, M. 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. doi:10.1243/09544062JMES1917

CrossRef Full Text | Google Scholar

Venkiteswaran, V. K., Sikorski, J., and Misra, S. (2019). Shape and Contact Force Estimation of Continuum Manipulators Using Pseudo Rigid Body Models. Mech. Mach. Theor. 139, 34–45. doi:10.1016/j.mechmachtheory.2019.04.008

CrossRef Full Text | Google Scholar

Wan, E. A., and Van Der Merwe, R. (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. doi:10.1109/ASSPCC.2000.882463

CrossRef Full Text | Google Scholar

Wu, J., Li, B., Li, C., Qiu, Y., Zhang, Z., Lyu, W., 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. doi:10.1117/12.2590684

CrossRef Full Text | Google Scholar

Xu, K., and Simaan, N. (2010). Intrinsic Wrench Estimation and its Performance index for Multisegment Continuum Robots. IEEE Trans. Robot. 26, 555–561. doi:10.1109/TRO.2010.2046924

CrossRef Full Text | Google Scholar

Yasin, R., and Simaan, N. (2020). Joint-level Force Sensing for Indirect Hybrid Force/position Control of Continuum Robots with Friction. Int. J. Rob. Res. doi:10.1177/0278364920979721

CrossRef Full Text | Google Scholar

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.

Edited by:

Chaoyang Shi, Tianjin University, China

Reviewed by:

JianXiong Hao, Tianjin University, China
Xiaoyang Li, Tianjin University, China

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