A Proposal of Haptic Technology to be Used in Medical Simulation

Abstract

For medical training aims, tele-operation systems have inspired virtual reality systems. Since force sensors placed on the robotic arms provide interaction force information that is transmitted to the human operator, such force produces a tactile sensation that allows feeling some remote or virtual environment properties. However, in the last two decades, researchers have focused on visually simulating the virtual environments present in a surgical environment. This implies that methods that cannot reproduce some characteristics of virtual surfaces, such as the case of penetrable objects, generate the force response. To solve this problem, we study a virtual reality system with haptic feedback using a tele-operation approach. By defining the operator-manipulated interface as the master robot and the virtual environment as the slave robot, we have, by addressing the virtual environment as a restricted motion problem, the force response. Therefore, we implement a control algorithm, based on a tele-operation system, to feedback the corresponding force to the operator. We achieve this through the design of a virtual environment using the dynamic model of the robot in contact with holonomic and non-holonomic constraints. In addition, according to the medical training simulator, before contact, there is always a free movement stage.

1. Introduction

Teleoperation and virtual reality systems are intrinsically related since they make a human operator interact with the environments without being in physical contact with them. In the first one, these environments are real, while in the second one, we generate them in a computer simulation. However, both types of systems must make the operator perceive, as realistically as possible, the characteristics of the remote or virtual environment. Some variables used to reproduce these characteristics are position and force, which provide visual and haptic feedback, respectively. The biggest difference between them is the procedure of generating the information received by the operator. In the teleoperation case, both signals exist physically and are transmitted via a control algorithm. The algorithm receives both signals from sensors. Both signals do not exist in virtual reality and we must generate computationally them.

Stimulating the senses of sight and touch, as precisely as possible, is essential during the interaction process, since they are the principal channels with which the operator perceives the world around them. For teleoperation, in the visual case, the communication comes directly from the environment or, if the operator is in a remote place, using a camera and a computer screen. The virtual reality system generates the environment through a digital simulation, and the operator receives the visual information through a screen. We need additional devices since the tactile issue is more complex. Those devices must be capable of transmitting the generated forces in the environment. Such a process implies including haptic robots in the systems because of their capability to generate forces and torques that the human operator can perceive in a tactile way. We need for teleoperation systems two physical robots, while in a virtual reality system we only need one robot and the virtual environment as the other.

The medical area has actively seized on both teleoperation and virtual reality systems. In the first case, a specialist can perform surgery procedures over long distances, eliminating the need for the physical presence of either the physicians or patients in the same location [1]. Practically, the specialist can examine or operate on the patient at a different geographic location without having to travel.

In Figure 1, we showed the emblematic Da Vinci Surgical Robotic System. Operating based on a master-slave control concept. The system provides the medical expert with a realistic operating environment that includes a high-quality stereo visualization and a human-machine interface that directly transfers the doctor’s hand gestures to the instrument inside the patient [2].

Figure 1.

Da Vinci surgical robotic system.

In the second case, someone has widely used virtual reality systems for medical training simulation. With the development of computer graphics, nowadays practically any surgical procedure can be visually simulated. Minimally invasive surgery has been the most beneficial area. It has implemented virtual environments in laparoscopy, neurosurgery, and urology, to name a few [3]. As an example, in Figure 2 we showed the simulator of Transurethral Resection of the Prostate TURP Mentor. Clinicians can improve their skills without endangering their patients to avoid Live practice with humans by using this system. Tactile feedback inclusion is an important factor in the improvement of such skills. It must synchronize tactile feedback with both the virtual reality simulation and the operator’s movements.

Figure 2.

URP Mentor (simulator).

The challenge for researchers in graphic computing and control systems is to design mathematical tools that fit with the object’s physical characteristics to be simulated within the virtual environment. Regarding the visual feedback, the position where such objects are located is essential for the operator to perceive his movements in the Cartesian virtual space. About the force, reproducing rigidity and softness takes special interest when the virtual environment includes penetrable and nonpenetrable objects. Here, the complexity of the mathematical tools increases because their physical laws are not always easy to simulate on a computer. For this reason, it must establish an interchange between visual and haptic realism because of the finite capacity of digital processing.

1.1 State-of-the-art

Teleoperation and virtual reality applications have developed in areas as different as automotive and video games. However, an important application is for medical surgery [4], where the operator on the master side needs to be sure of the force he feels. Such a force, generated in opposition to his movements, must be ideally the same as that the robot applies over the patient at the slave side. It has widely used virtual reality systems in minimally invasive surgical simulation, where the operator should feel the same forces that it would feel in a real procedure, [5].

Goertz presented the first master-slave teleoperation system [6]. It was used to handle toxic waste using two coupled manipulators. Subsequently, using electrical signals, he included a rudimentary force generation device. Since then, we have used this system in areas such as micro and nano-manipulation [7] underwater exploration [8], and tele-surgery [9]. Haptic feedback takes special relevance in this last area since it is crucial for the surgeon to receive an accurate force response. The most effective way to do so is to have force sensors at both master and slave sides, and a control algorithm that guarantees an accurate tracking between the contact forces present in the remote environment and those sent to the operator [10].

The idea of force feedback on virtual reality systems began with the fundamental work of Sutherland [11]. He established that the interaction between the human operator and the virtual environment should not only be visual but also tactile. It was not until the 1990s that he adapted the Goertz device to provide force feedback during virtual molecular coupling [12]. Since then, the use of manipulators in virtual reality applications has spread to CAD/CAM assembly [13], aerospace maintenance [14], and especially in medical training through simulators [15] where, unlike systems of master-slave teleoperation, neither the environments nor the contact forces exist. It must transmit the actual forces to the operator with precision. The quality of this transmission depends on the characteristics of the haptic interface and the corresponding force control algorithm [3].

Articulated robots play a major role in medical training simulation systems with force feedback since this kind of electromechanical device can measure spatial position and generate torques. There has been a large effort to design robot haptic interfaces such as the widely used Phantom serial robot [16], the Delta Haptic parallel robot [17], and the combination of passive elements such as brakes and springs with motors [18]. Such robots are examples of impedance types devices, i.e. they read position and control force in response. There are another robots that read forces and control motion, called admittance type devices. The difference between using one or another type relies on the characteristics of the virtual environment (e.g., stiffness, inertia, damping, friction).

Along with haptic interfaces, there has been an intense development of graphical simulation tools capable of reproducing a wide range of virtual environments. The principal aim is for the operator to perceive, as realistic as possible, objects with a high quality of detail. The applications developed include microscopic exploration [19], aviation [20], and clinical neuropsychology [21], among many others. With medical training simulation, a correct synergy between visual and tactile feedback is essential to heighten the skills of medicine students. However, to increase immersion and consequently the realism of the virtual reality displays, it is necessary to model environments that combine haptic and graphics to the same complexity [22]. This is not always possible since it limited the computational processing and it cannot execute the applications in real-time.

Salisbury et al. [23] presented the basic architecture for a virtual reality application with visual and haptic feedback. They established that the force rendering algorithms must be geometry-dependent. This is a disadvantage in medical training simulation since the virtual objects to be reproduced include bones and organs with irregularities or indentations. There are cases in which the interaction occurs not only on the surfaces of the object, but we must also calculate the penetration forces, as we do in surgery simulators. The alternative is to design algorithms based on physical laws that involve the dynamic and movement of the objects when the operator interacts with them [24]. The perfect scenario would be to render forces by combining physical approaches with the most sophisticated haptic interfaces. However, as mentioned before, doing that is computationally more expensive, not to mention the high costs it would entail.

The more realistic the force transmitted to the operator, the higher the quality of the method used. The factors mentioned above cause a series of compensations between haptic and visual realism, real-time execution, and system costs. Such offsets allowed for establishing two principal methods for rendering forces from virtual environments [25]. The first one is the penalty method, which is widely used for its simplicity since it required a penetration measure starting at the contact point with the virtual object. The second approach is the imposed motion method, where it is considering the contact as a bilateral constraint and it is calculating the contact force response using Lagrange multipliers. Computer graphics and haptic applications use both methods. Ruspini et al. established the differences between using both methods [22].

It is to assume, according to the penalty method, virtual objects are formed geometrically or algebraically by defined primitives such as lines, planes, spheres, and cylinders [26]. Therefore, the force rendering depends on an implicit equation and a contact point given by a collision detection algorithm. In 1991, Sclaroff and Pentland proposed a generalization using the implicit function representation method to allow collision detection for common 3D shapes [27]. The method is to replace the polygon and spline with their previous versions. By using this technique, it is possible to use local gradients in the normal direction of the virtual surface [28]. In this sense, the concept of impedance takes special relevance since it addressed the force rendering problem as an energy exchange phenomenon, allowing to study of the stability of the haptic system [29].

Using the imposed motion method implies considering the contacts with the virtual object as bilateral constraints [25]. Therefore, haptic systems are a master-slave scheme where the rigid object, the virtual representation of the haptic interface, is the master. This entails employing Lagrange multipliers to compute the magnitude of the contact force. In 1994, Bayo and Avello designed an algorithm by considering the dynamics of a multi-body system as a constraint [30]. The advantage of this approach is that the force response can not only be in terms of a single contact point, but it can also be in terms of the dynamic characteristics of the surfaces [31]. However, the hardware requirement for haptic rendering in 3D is for the haptic interface to have at least 3 degrees of freedom.

The methods mentioned above have been the cornerstone for virtual forces generation both in graphic computing and haptic systems. One of the principal requirements in such areas is that the systems be capable of reproducing the forces that would be present during contact with rigid and soft objects. This is especially important in the development of simulators used in medical training, where the tactile sensations caused by contact between a virtual tool and bones or organs must reproduce [32]. Depending on the goals of the design, we must make a compromise between haptic and visual realism.

2. Preliminaries

In order to implement a virtual reality application, it is important to combine and match both visual and haptic feedback in real-time. Next, we introduce the fundamentals of virtual surface representation, which allow us to show how works the force rendering method, proposed in [33]. The mathematical nature of holonomic and non-holonomic constraints is presented by introducing some basic concepts of differential topology and the mathematical model of the haptic system. For this, a teloperation scheme is adapted considering the slave robot as virtual. It simulated its dynamic model within the virtual environment. Said robot is in contact with holonomic or non-holonomic virtual constraints, whose mathematical representation is also presented. The process for rendering the forces to be transmitted to the operator is detailed below, as well as some observations on the validity of the approach, especially regarding the differences when using both types of constraints.

2.1 Virtual surfaces representation

It is important to distinguish between the different representations of virtual surfaces, both from the point of view of computer graphics and from the point of view of haptic systems. Because the physics of light (with visual representation) differs from the physics of mechanical interactions. It is important to consider that although the graphical and haptic simulation can share the coding of certain properties, such as the shape, they must differ in other aspects, such as models, mathematical techniques, and implementation [32]. It is important to note that haptic rendering in practice avoids the many complex renderings developed by the graphics computing community over the years. However, in this section, the basics of such representations are given to introduce the fundamentals of the proposed haptic representation method, which central idea is to relate holonomic and non-holonomic constraints, which mathematically have a kinematic basis, with the rigid and soft tissues of the body. Surface dynamic complexity is avoided, and instead, from a purely haptic approach, it is classified as nonpenetrable and penetrable. By using the approach of a manipulator in constrained motion, it includes the dynamic of the virtual robot to exemplify that it is possible to model the virtual tool, in a more complex way than that of a single point probe used commonly in computer graphics [34].

2.1.1 Nonpenetrable virtual surfaces

In order to represent forces coming from rigid surfaces, it suffices to define algebraically an implicit equation in joint coordinates (task space), or at least two of its geometric characteristics such as normal and tangential vectors, or distance and angle relationships between points, lines, and planes [35]. This simplifies the graphical and haptic implementation since it is possible to define the surface as the zero set of a function $f$ valued at $\mathbb{R}$ as $S_f=\left\{\mathit{x}\in {\mathbb{R}}^3|f\left(\mathit{x}\right)=0\right\}$ [36]. Based on this approach, lines, points, and mainly zero-width polygons, assembled by vertices, form rigid virtual environments, [31]. For example, in Figure 3, a human skull using a triangle mesh with Phys X${}^{\circledR }$ by NVIDIA graphics engine is modeled.

Figure 3.

Convex decomposition by Phys X${}^{\circledR }$. (a) Triangle mesh collision. (b) Convex decomposition.

In this case, a PxGeometry Class by NVIDIA, part of a geometry class, from the common base class, is used. Each geometry class defines a volume or surface with a fixed position and orientation. We can implement it in many geometry types, as simple as spheres, capsules, boxes and planes, and others more complex as convex meshes, triangle meshes, and height fields [37]. Besides, the methods to build a complex object as a human skull and its properties include triangle mesh collision and convex decomposition, Figure 3. The haptic interaction for virtual objects formed by such geometries and methods comes from a normal vector over the geometry, as we can see in Figure 4, where the haptic interaction with a sphere occurs using a collision detection algorithm of a single point.

Figure 4.

Haptic interaction with a sphere. (a) Applied external forces. (b) Collision detection. (c) Collision response. (a) Original mesh. (b) Delaunay triangulation. (c) Mass-spring like system.

2.1.2 Penetrable virtual surfaces

The case of deformable surfaces is more complex since, from a biomedical approach, a physically realistic simulation must consider all the nonlinearities of material deformation (e.g. stress, strain, elasticity, and viscoelasticity). One strategy is to combine a finite element discretization of the geometry together with a finite difference discretization of time and an updated Lagrangian iterative scheme [38]. Another very used representations of deformable surfaces in computer graphics are the particle-based model. Their location, speed, acceleration, mass, and any other parameter needed for an application describe particles and they develop according to Newtonian mechanical laws [36]. For example, in Figure 5, we model a deformable tissue using the graphic engine FLEX${}^{\circledR }$ by NVIDIA where, from a polygon-based mesh, a particle system is got through a Delaunay triangulation.

Figure 5.

Particle model by FLEX (soft tissue). (a) Original mesh. (b) Delaunay triangulation. (c) Mass-spring like system.

Using particle-based models allows to visually reproduce more complex processes such as cutting and slitting, processes very common in medical training. However, no matter how complicated the underlying model is, the force response because of deformation is a function of deflection only, that is, we define deformation as the displacement of the initial contact point between an instrument and a body. Deformable, even if it is a very large deformation, [32]. This characteristic allows us to study the contact from a kinematic perspective, it moved once the manner tool inside the soft object.

More sophisticated graphic tools such as SOFA (Simulation Open Framework Architecture) address the description of the object, typically by using three models: an internal model with independent degrees of freedom, the mass, and the constitutive laws, a collision model with contact geometry, and a visual model with detailed geometry and rendering parameters [39]. During runtime, it synchronizes the models using a generic mechanism called mapping, which is used to propagate forces and displacements and to enforce coherence between them. Normally, the internal model acts as the master system, imposing its displacements on the slave systems (the visual and collision model).

Let $f$ be the function used to map the positions $x_{\mathrm{m}}$ of a master model to the position $x_{\mathrm{s}}$ of a slave model

$${\mathit{x}}_{\mathrm{s}}=f\left({\mathit{x}}_{\mathrm{m}}\right).$$E1

We map the velocities similarly to

$${\dot{\mathit{x}}}_{\mathrm{s}}=\mathit{J}\left({\mathit{x}}_{\mathrm{m}}\right){\dot{\mathit{x}}}_{\mathrm{m}}$$E2

where the Jacobian matrix $\mathit{J}\left({\mathit{x}}_{\mathrm{m}}\right)$ encodes the linear relationship between the velocities of the master and slave systems. In linear assignments, the operators $f$ and $\mathit{J}\left({\mathit{x}}_{\mathrm{m}}\right)$ are the same, otherwise $f$ is not linear regarding ${\mathit{x}}_{\mathrm{m}}$, and we cannot write it as a matrix.

Given forces ${\mathit{\lambda }}_{\mathrm{s}}$ applied to a slave model, the mapping computes the equivalent forces ${\mathit{\lambda }}_{\mathrm{m}}$ applied to its master. Since equivalent forces must have the same energy, [39], the following relation holds

$${\dot{\mathit{x}}}_{\mathrm{m}}^{\mathrm{T}}{\mathit{\lambda }}_{\mathrm{m}}={\dot{\mathit{x}}}_{\mathrm{s}}^{\mathrm{T}}{\mathit{\lambda }}_{\mathrm{s}}.$$E3

The kinematic relation (2) allows to rewrite the Eq. (3) as

$${\dot{\mathit{x}}}_{\mathrm{m}}^{\mathrm{T}}{\mathit{\lambda }}_{\mathrm{m}}={\mathit{J}}^{\mathrm{T}}\left({\mathit{x}}_{\mathrm{m}}\right){\dot{\mathit{x}}}_{\mathrm{m}}^{\mathrm{T}}{\mathit{\lambda }}_{\mathrm{s}}.$$E4

Since Eq. (4) holds for all velocities ${\dot{\mathit{x}}}_{\mathrm{m}}$, the principle of virtual work allows us to simplify it to get

$${\mathit{\lambda }}_{\mathrm{m}}={\mathit{J}}^{\mathrm{T}}\left({\mathit{x}}_{\mathrm{m}}\right){\mathit{\lambda }}_{\mathrm{s}}.$$E5

The kinematic mappings (1), (2), and (5) allow to compute displacements and to apply forces. They are also used to connect generalized coordinates, such as joint angles, to task space geometries.

2.2 Geometry of a constrained sub-manifold

From a robot control approach, we can consider a tool in contact with an object as a robot in constrained motion. The constraints of this system will be well defined if we associate them with physically realizable forces. This occurs, for example, with an industrial robot in contact with a proper surface (a real one) like a car bonnet in a painting or welding task. But with virtual environments, where surfaces do not exist, there are no physical constraint forces associated to them. Thus, the constraints are not well defined, and they are called virtual constraints, [40]. It mathematically addressed the nonpenetrable and penetrable virtual surfaces as virtual constraints. In this sense, it is important to introduce the geometric properties of such constraints in order to define them as either holonomic or non-holonomic.

Let $Q$ be the $n$-dimensional smooth manifold configuration space of an unconstrained manipulator and $\mathit{q}\in {\mathbb{R}}^n$ its local generalized coordinates. The tangent space to $Q$ at $\mathit{q}$, denoted ${\mathit{\tau }}_{\mathit{q}}Q$ comprises all generalized velocity vectors $\dot{\mathit{q}}\in {\mathbb{R}}^n$ of the system.

Definition 2.1. A geometric constraint on $Q$ is a relation of the form

$${\mathit{h}}_i\left(\mathit{q}\right)=\mathbf{0}i=1,\dots ,k<n,$$E6

where ${\mathit{h}}_i:Q\to \mathbb{R}$ limits the admissible motions of the system to a $\left(n-k\right)$-dimensional smooth sub-manifold of $Q$■.

Constraints that involve not only the generalized coordinates but also their first derivatives in the form

$${\mathit{a}}_i\left(\mathit{q},\dot{\mathit{q}}\right)=\mathbf{0}i=1,\dots ,k<n,$$E7

with ${\mathit{a}}_i\left(\mathit{q},\dot{\mathit{q}}\right)\in {\mathit{\tau }}_{\mathit{q}}Q$, are called kinematic constraints. These limit the allowable movements of the manipulator to a $\left(n-k\right)$-dimensional smooth sub-manifold of $Q$ by restricting the set of generalized velocities that can be achieved in a configuration.

Definition 2.2. A Pfaffian constraint on $Q$ is a set of $k$ kinematic constraints, which are linear in velocity in the following form

$${\mathit{a}}_i^{\mathrm{T}}\left(\mathit{q},\dot{\mathit{q}}\right)=\mathbf{0}i=1,\dots ,k<n,$$E8

where ${\mathit{a}}_i\left(\mathit{q},\dot{\mathit{q}}\right):Q\to {\mathrm{R}}^n$ is assumed to be smooth and linearly independent.■

A kinematic constraint can be integrable, there are $k$ real-valued functions ${\mathit{h}}_i\left(\mathit{q}\right)$ such that

$$\frac{\partial }{\partial \mathit{q}}\left({\mathit{h}}_i\left(\mathit{q}\right)\right)={\mathit{a}}_i^{\mathrm{T}}\left(\mathit{q}\right)i=1,\dots ,k<n.$$E9

Here, the kinematic constraints are, in fact, geometric constraints. Pfaffian constraints set of ${\mathit{a}}_i\left(\mathit{q}\right)$ is known as holonomic if it is integrable, the system has a geometric limitation. For example, by considering a set of holonomic constraints characterized by

$${\mathit{\phi }}_i\left(\mathit{q}\right)=\mathbf{0}i=1,\dots ,k<n.$$E10

From Eq. (9), we get

$$\frac{\partial }{\partial \mathit{q}}\left({\mathit{\phi }}_i\left(\mathit{q}\right)\right)={\mathit{J}}_{\phi }^{\mathrm{T}}\left(\mathit{q}\right)$$E11

where ${\mathit{J}}_{\phi }\left(\mathit{q}\right)\in {\mathbb{R}}^{k\times n}$ is the Jacobian of the holonomic constraint. Therefore, holonomic constraints are characterized by equivalent equations in terms of position variables, we can get the position equations by integrating them, if velocity equations initially described the constraints, [41].

Property 2.1. Given $n$ generalized coordinates $\mathit{q}$ in a sub-manifold $Q$ and $k$ holonomic constraints, the space tangent to $Q$ in a configuration can be described by adequately defining $\left(n-k\right)$ new generalized coordinates of the restricted sub-manifold that characterize the real degrees of freedom of the system [42].■

The set of Pfaffian constraints ${\mathit{a}}_i\left(\mathit{q}\right)$ is called non-holonomic if it is nonintegrable, the system has a kinematic limitation. Assuming again that the vectors ${\mathit{a}}_i:Q\to {\mathbb{R}}^n$ are smooth and linearly independent, the non-holonomic constraints can be expressed as

$$\mathit{A}\left(\mathit{q}\right)\dot{\mathit{q}}=\mathbf{0}$$E12

where $\mathit{A}\left(\mathit{q}\right)\in {\mathbb{R}}^{k\times n}$ is the Pfaffian matrix of non-holonomic constraints and which image space produces forces to ensure that the system does not move in those directions. The presence of these constraints limits the system mobility completely differently if compared to holonomic ones, even if its generalized velocities at each point are constrained to a $\left(n-k\right)$ dimensional sub-manifold space, it is still possible to reach any configuration in $Q$.

Property 2.2. Given $n$ generalized coordinates $\mathit{q}$ in a sub-manifold $Q$ and $k$ non-holonomic constraints, the space tangent to $Q$ in a configuration has $\left(n-k\right)$ degrees of freedom but the number of generalized coordinates cannot be reduced, [43].■

Remark 2.1. We assume that velocity-level gives non-holonomic constraints Eq. (12), and position-levelEq. (10) describes holonomic constraints. In practical problems, it may describe both types of constraints as velocity-level equations.■

2.2.1 Integrability of the constraints

A vector field $\mathit{g}:{\mathbb{R}}^n\to {\mathit{\tau }}_{\mathit{q}}{\mathbb{R}}^n$ is a smooth mapping assigning to each point $\mathit{q}\in {\mathbb{R}}^n$ a tangent vector $\mathit{g}\left(\mathit{q}\right)\in {\mathit{\tau }}_{\mathit{q}}{\mathbb{R}}^n$. In local coordinates, we can represent $\mathit{q}$ as a column vector whose elements depend on $\mathit{q}$ as

$$\mathit{g}\left(\mathit{q}\right)=\left[\begin{array}{c}{\mathit{g}}_1\left(\mathit{q}\right)\\ ⋮\\ {\mathit{g}}_n\left(\mathit{q}\right)\end{array}\right],$$E13

where $\mathit{g}$ is smooth if each ${\mathit{g}}_i\left(\mathit{q}\right)$ is smooth.

Given ${\mathit{g}}_1$ and ${\mathit{g}}_2$, we define the Lie bracket of these vectors fields as

$$\left[{\mathit{g}}_1,{\mathit{g}}_2\right]=\frac{\partial }{\partial \mathit{q}}\left({\mathit{g}}_2\right){\mathit{g}}_1-\frac{\partial }{\partial \mathit{q}}\left({\mathit{g}}_1\right){\mathit{g}}_2$$E14

where $\left[{\mathit{g}}_1,{\mathit{g}}_2\right]$ is a new vector field.

A distribution assigns a subspace of the tangent space to each point in ${\mathbb{R}}^n$ smoothly. A special case is a distribution defined by a set of smooth vector fields, ${\mathit{g}}_1,\dots ,{\mathit{g}}_m$. Here, we can define distribution as

$$\mathrm{\Delta }=\text{span}\left\{{\mathit{g}}_1,\dots ,{\mathit{g}}_m\right\},$$E15

where the span over the set of smooth real-valued functions on ${\mathbb{R}}^n$ is taken. Evaluated at any point $\mathit{q}\in {\mathbb{R}}^n$, the distribution defines a linear subspace of the tangent space

$${\mathrm{\Delta }}_{\mathit{q}}=\text{span}\left\{{\mathit{g}}_1\left(\mathit{q}\right),\dots ,{\mathit{g}}_m\left(\mathit{q}\right)\right\}\subset {\mathit{\tau }}_{\mathit{q}}{\mathbb{R}}^n.$$E16

A distribution is involutive if it is closed under the Lie bracket, i.e.,

$$\left[{\mathit{g}}_i,{\mathit{g}}_j\right]\in \mathrm{\Delta },\forall {\mathit{g}}_i,{\mathit{g}}_j\in \mathrm{\Delta }.$$E17

We said distribution $\mathrm{\Delta }$ of dimension $k$ to be integrable if, for every point $\mathit{q}\in {\mathbb{R}}^n$, there are a set of smooth functions ${\mathit{h}}_i:{\mathbb{R}}^n\to \mathbb{R}$ for $i=1,\dots ,n-k$ such that the row vectors $\frac{\partial }{\partial \mathit{q}}\left({\mathit{h}}_i\right)$ are linearly independent at $\mathbf{q}$ and for every $\mathbf{g}\in \mathrm{\Delta }$.

$$\frac{\partial }{\partial \mathit{q}}\left({\mathit{h}}_i\right)\mathit{g}\left(\mathit{q}\right)=\mathbf{0}i=1,\dots ,n-k.$$E18

The hypersurfaces defined by the level sets $\left\{{\mathit{h}}_1\left(\mathit{q}\right)=c_1,\dots ,{\mathit{h}}_{n-k}\left(\mathit{q}\right)=c_{n-k}\right\}$ are called integral manifolds for the distribution $\mathrm{\Delta }$. Eq. (18) shows that $\mathrm{\Delta }$ coincides with the tangent space to its integral manifold at $\mathit{q}$. We relate integral manifolds to involutive distributions by the following so-called Frobenious theorem, [43].

Theorem 2.1. A distribution is integrable if and only if it is involutive.■

This theorem gives a necessary and sufficient condition for the complete integrability of a distribution. Thus, if $\mathrm{\Delta }$ is a $k$-dimensional involutive distribution, then locally there are $n-k$ functions ${\mathit{h}}_i:{\mathbb{R}}^n\to \mathbb{R}$ such the level surfaces $\mathit{h}=\left({\mathit{h}}_1,\dots ,{\mathit{h}}_{n-k}\right)$ give that integral manifolds of $\mathrm{\Delta }$. The result mentioned above gives conditions for the integrability of a set of kinematic constraints in the following proposition.

Proposition 2.1. [44] The set of $k$ Pfaffian constraints, described in (8), is holonomic if and only if $\mathrm{\Delta }$ its an involutive distribution.■

So it is possible to establish when a Pfaffian constraint is non-holonomic by checking if its distribution is not involutive.

2.3 Haptic system overview

A common practice in the computer graphics community is to associate the position and orientation of virtual tools directly with that of the haptic interface. However, this assumption is because some real tools have negligible dynamics, such as a scalpel. From a teleoperation approach, we assumed that even the simplest tool has some dynamic properties to consider in the virtual environment. This section presents a description of this proposal, both mathematical and intuitive.

In order to describe the operation of the haptic system, two independent sets of task space coordinates are considered as shown in Figure 6.

Figure 6.

Haptic system.

The operator manipulates the haptic interface, i.e., the master robot in the real environment and we denote whose Cartesian coordinates as ${\mathit{x}}_{\mathrm{m}}\in \mathit{SE}\left(3\right)$, where ${\mathit{p}}_{\mathrm{m}}\in {\mathrm{R}}^3$ is the end-effector position and ${\mathit{R}}_{\mathrm{m}}\in \mathit{SO}\left(3\right)$ its orientation. The virtual tool must respond to the movements of such interface in the virtual environment with Cartesian coordinates ${\mathit{x}}_{\mathrm{v}}\in \mathit{SE}\left(3\right)$, where ${\mathit{p}}_{\mathrm{v}}\in {\mathbb{R}}^3$ is the virtual tool position and ${\mathit{R}}_{\mathrm{v}}\in \mathit{SO}\left(3\right)$ its orientation. In a teleoperation context, the position of the master robot acts as a reference for the virtual tool and it is visually projected on the screen through the virtual avatar of the system. The operator moves the virtual tool freely until the collision detection algorithm shows that contact with the virtual surface is taking place. Until now, the master robot exerts a force that is measured by a force sensor that serves as a reference for the virtual robot and must apply to that surface. By closing the feedback loop, the control algorithm produces a tactile sensation for the operator. Ideally, both visual and haptic feedback must coincide, allowing the operator to have a visual reference to the virtual tool and the feeling of the dynamic changes of its contact with the virtual surface.

In a similar approach, Faurlin et al. clarify in [44] that the virtual environment can be represented by a set of generalized coordinates ${\mathit{q}}_{\mathrm{v}}\in {\mathbb{R}}^3$, which are related to the task-space coordinates of the master robot by a nonlinear kinematic equation

$${\mathit{x}}_{\mathrm{m}}=f\left({\mathit{q}}_{\mathrm{v}}\right)$$E19

which is a mapping between the real and virtual environment, similar to that of Eq. (1) but where the former acts as a master model.

The set of coordinates ${\mathit{q}}_{\mathrm{v}}$ allows the dynamic model of the virtual tool to be described in terms of Euler-Lagrange equations of motion. A set of holonomic or non-holonomic constraints that represent the virtual surface can be embedded into the kinematic mapping (19), relating independent master robot task-space coordinates and dependent virtual robot task-space coordinates. The virtual tool moves according to the physic simulation propagated in the virtual environment coordinates and always satisfies such constraints.

2.3.1 Dynamic model and properties

Consider a real master (m) and a virtual slave (v) robot system composed of two manipulators, each of them with $n$ degrees of freedom but not necessarily with the same kinematic configuration. Each robot spans a $k$-dimensional task space and, based on master/virtual devices, can be scaled to meet the desired virtual application. The master dynamics is given by

$${\mathit{H}}_{\mathrm{m}}\left({\mathit{q}}_{\mathrm{m}}\right){\ddot{\mathit{q}}}_{\mathrm{m}}+{\mathit{C}}_{\mathrm{m}}\left({\mathit{q}}_{\mathrm{m}},{\dot{\mathit{q}}}_{\mathrm{m}}\right){\dot{\mathit{q}}}_{\mathrm{m}}+{\mathit{D}}_{\mathrm{m}}{\dot{\mathit{q}}}_{\mathrm{m}}+{\mathit{g}}_{\mathrm{m}}\left({\mathit{q}}_{\mathrm{m}}\right)={\mathit{\tau }}_{\mathrm{m}}-{\mathit{\tau }}_{\mathrm{h}}$$E20

while the virtual slave dynamics is modeled by

$${\mathit{H}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right){\ddot{\mathit{q}}}_{\mathrm{v}}+{\mathit{C}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}},{\dot{\mathit{q}}}_{\mathrm{v}}\right){\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{D}}_{\mathrm{v}}{\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{g}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)={\mathit{\tau }}_{\mathrm{v}}+{\mathit{\tau }}_{\mathrm{s}}$$E21

where the subscripts m and v denote the real master and the virtual slave manipulators, respectively. For $i=m,v,{\mathit{q}}_i\in {\mathbb{R}}^n$ is the vector of generalized coordinates, ${\mathit{H}}_i\left({\mathit{q}}_i\right)\in {\mathbb{R}}^{n\times n}$ is the inertia matrix, ${\mathit{C}}_i\left({\mathit{q}}_i,{\dot{\mathit{q}}}_i\right){\dot{\mathit{q}}}_i\in {\mathbb{R}}^{n\times n}$ is the vector of Coriolis and centripetal forces, ${\mathit{D}}_i{\dot{\mathit{q}}}_i\in {\mathbb{R}}^{n\times n}$ is a diagonal matrix of viscous friction coefficients, ${\mathit{g}}_i\left({\mathit{q}}_i\right)\in {\mathbb{R}}^n$ is the vector of gravitational torques, ${\mathit{\tau }}_i\in {\mathbb{R}}^n$ is the vector of generalized inputs, ${\mathit{\tau }}_{\mathrm{h}}\in {\mathbb{R}}^n$ is the real torque applied by the human operator on the master side and ${\mathit{\tau }}_{\mathrm{s}}\in {\mathbb{R}}^n$ is the virtual torque generated because of the contact with the virtual constraint [45].

Property 2.3. With a proper definition of the robot parameters, it is possible to express the robot dynamics as

$${\mathit{H}}_i\left({\mathit{q}}_i\right){\ddot{\mathit{q}}}_i+{\mathit{C}}_i\left({\mathit{q}}_i,{\dot{\mathit{q}}}_i\right){\dot{\mathit{q}}}_i+{\mathit{D}}_i{\dot{\mathit{q}}}_i+{\mathit{g}}_i\left({\mathit{q}}_i\right)={\mathit{Y}}_i\left({\mathit{q}}_i,{\dot{\mathit{q}}}_i,{\ddot{\mathit{q}}}_i\right){\mathit{\theta }}_i$$E22

where ${\mathit{Y}}_i\left({\mathit{q}}_i,{\dot{\mathit{q}}}_i,{\ddot{\mathit{q}}}_i\right)\in {\mathbb{R}}^{n_i\times l}$ is the regressor and ${\mathit{\theta }}_i\in {\mathbb{R}}^l$ is a constant vector of parameters.■

Assumption 2.1. The master and the virtual slave robots share the same geometric structure, but they do not necessarily have the same parameters of the dynamic model, that is, the matrices and vectors of the models described in (20) and (21) do not have to be the same.■

External torques are acting in both robots, either the real torque ${\mathit{\tau }}_{\mathrm{h}}$ applied by the human on the master side or the virtual torque ${\mathit{\tau }}_{\mathrm{s}}$ generated because of the contact between the virtual robot and the virtual surface. We can define the torque applied by the human operator as

$${\mathit{\tau }}_{\mathrm{h}}={\mathit{J}}_{\mathrm{m}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{m}}\right){\mathit{F}}_{\mathrm{h}}$$E23

where ${\mathit{F}}_{\mathrm{h}}\in {\mathbb{R}}^3$ is the force applied by the operator in the task-space coordinates and ${\mathit{J}}_{\mathrm{m}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{m}}\right)\in {\mathbb{R}}^{3\times n}$ the geometric Jacobian of the master manipulator. In the same way, the torque applied on the virtual surface can be expressed as

$${\mathit{\tau }}_{\mathrm{s}}={\mathit{J}}_{\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right){\mathit{F}}_{\mathrm{s}}$$E24

where ${\mathit{F}}_{\mathrm{v}}\in {\mathbb{R}}^3$ is the force applied on such surface in task-space coordinates.

2.3.2 Virtual holonomic constraints

Whit of holonomic constraints we assume that, in virtual task space coordinates, the virtual robot is subject to $k$ virtual holonomic constraints characterized by

$${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)=\mathbf{0}$$E25

where a suitable normalization is done for the gradient of this constraint,

$${\mathit{J}}_{\phi x_{\mathrm{v}}}\left({\mathit{x}}_{\mathrm{v}}\right)=\nabla {\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)\in {\mathbb{R}}^{k\times n},$$E26

to be unitary.

The representation of constraint (25) in generalized virtual coordinates is

$${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)=\mathbf{0}$$E27

where ${\mathit{q}}_{\mathrm{v}}\in {\mathbb{R}}^n$ is the vector of the virtual robot end-effector joint coordinates. The gradient of the constraint (27) is

$${\mathit{J}}_{\phi \mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)=\nabla {\mathit{\phi }}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)\in {\mathbb{R}}^{k\times n}.$$E28

These two gradients are related by

$${\mathit{J}}_{{\mathit{\phi x}}_{\mathrm{v}}}\left({\mathit{q}}_{\mathrm{v}}\right)={\mathit{J}}_{{\mathit{\phi x}}_{\mathrm{v}}}\left({\mathit{x}}_{\mathrm{v}}\right){\mathit{J}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)$$E29

where ${\mathit{J}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)\in {\mathbb{R}}^k$ is the geometric Jacobian of the virtual manipulator. Hence, the torque because of the contact with the virtual surface in (1) can be defined as

$${\mathit{\tau }}_{\mathrm{s}}={\mathit{J}}_{\mathit{\phi }\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right){\mathit{\lambda }}_{\mathrm{v}}$$E30

where ${\mathit{\lambda }}_{\mathrm{v}}\in {\mathbb{R}}^k$ is a vector of Lagrange multipliers that represents the virtual force applied over the surface. Then, it is possible to rewrite the whole Eq. (21) as

$${\mathit{H}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right){\ddot{\mathit{q}}}_{\mathrm{v}}+{\mathit{C}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}},{\dot{\mathit{q}}}_{\mathrm{v}}\right){\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{D}}_{\mathrm{v}}{\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{g}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)={\mathit{\tau }}_{\mathrm{v}}+{\mathit{J}}_{\mathit{\phi }\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right){\mathit{\lambda }}_{\mathrm{v}}.$$E31

According to Property 2.1, the virtual holonomic constraints (2.26) reduce the number of degrees of freedom of the virtual robot and the dimension of its configuration space to a $n-k$–dimensional sub-manifold, [43].

2.3.3 Virtual non-holonomic constraints

With non-holonomic constraints, something well known is that we cannot express them as a function of only the generalized coordinates as in (25) or (27). Instead, they are commonly expressed as Pfaffian constraints. In the present case, these constraints are written more intuitively in terms of the virtual end-effector velocities ${\mathit{v}}_{\mathrm{v}}={\left[{\dot{\mathit{p}}}_{\mathrm{v}}{\mathit{\omega }}_{\mathrm{v}}\right]}^{\mathrm{T}}$ as

$${\mathit{A}}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right){\mathit{v}}_{\mathrm{v}}=\mathbf{0},$$E32

where ${\dot{\mathit{p}}}_{\mathrm{v}}\in {\mathbb{R}}^3$ and ${\mathit{\omega }}_{\mathrm{v}}\in {\mathbb{R}}^3$ are the linear and angular velocities of the virtual end-effector, respectively, and ${\mathit{A}}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)\in {\mathbb{R}}^{k\times n}$ is the corresponding Pfaffian constraint matrix. If the dynamic equations are defined in the virtual joint-space coordinates ${\mathit{q}}_{\mathrm{v}}$, these constraints are projected via Faurling et al. [44].

$${\mathit{A}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)={\mathit{A}}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right){\mathit{J}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)$$E33

Assuming that the virtual robot is subject to $k$ velocity-level equations of non-holonomic constraints characterized by

$${\mathit{A}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right){\dot{\mathit{q}}}_{\mathrm{v}}=\mathbf{0}$$E34

the torque because of the contact with the virtual environment in (21) can be expressed as

$${\mathit{\tau }}_{\mathrm{s}}={\mathit{A}}_{\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right){\mathit{\lambda }}_{\mathrm{v}}$$E35

where ${\mathit{\lambda }}_{\mathrm{v}}\in {\mathbb{R}}^k$ is the vector of Lagrange multipliers which determines the magnitude of the constraint forces over the virtual surface. Then, it is possible to rewrite Eq. (21) as

$${\mathit{H}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right){\ddot{\mathit{q}}}_{\mathrm{v}}+{\mathit{C}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}},{\dot{\mathit{q}}}_{\mathrm{v}}\right){\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{D}}_{\mathrm{v}}{\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{g}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)={\mathit{\tau }}_{\mathrm{v}}+{\mathit{A}}_{\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right){\mathit{\lambda }}_{\mathrm{v}}.$$E36

The non-holonomic constraints reduce the number of virtual robot available degrees of freedom to an $\left(n-k\right)$-dimensional sub-manifold, but they do not reduce the dimension of its configuration space [43, 46].

3. System implementation

In this section, the theoretical and practical aspects of implementing a virtual reality system with virtual restrictions are presented. The principal aspect concerned is the design of a controller capable to perform accurate haptic feedback that makes to feel the operator to be in contact with either a penetrable or nonpenetrable virtual surface. The method for visually reproducing the virtual tool in contact with the virtual objects is presented. It is important to note that this method avoids the complexity of the virtual environments currently implemented in the simulators used in medical training. However, the basic aspects addressed to show that the virtual constraints approach can be used practically, and eventually adapted to sophisticated graphic computing tools.

3.1 Virtual environment design

As mentioned in Section 2, the important aspect to get realistic haptic feedback from a surface embedded into a virtual environment comprises defining its geometry. In Figure 7, we show an idealized representation of a virtual point probe in contact with either an nonpenetrable or penetrable virtual surface.

Figure 7.

($\cdots$) Penetrable and (__) nonpenetrable surfaces.

In the first case, we assume that the contact arises at a single point over the surface from where the virtual probe cannot move forward, i.e., its velocity is equal to zero. Therefore, a normal force vector, which magnitude increases depending on the force applied by the operator, avoids motion. In robot control, if we are supposed to connect the probe to the end effector of a manipulator, according to Property 2.1, the number of degrees of freedom of the system in contact with the surface is reduced. Intuitively, that means that the virtual probe cannot move forward from where the contact arises, which can be at any point on the surface [47]. Actually, if the virtual object is built up by using a polygonal method, there will be a set of surfaces $\left({\mathit{\phi }}_0,{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_n\right)$ joined by vertices as shown in Figure 7.

The best way to find the place of the contact point (which belongs to a set of points defining each surface) is by establishing an implicit equation ${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)$ containing such a point. The set of points defining the virtual surface are

$${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)=\mathbf{0},$$E37

which coincides with the holonomic constraint of Eq. (25) in Cartesian coordinates or Eq. (27) in generalized coordinates. From those expressions, it can establish a collision detection algorithm by defining the following conditions:

• If ${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)>\mathbf{0}$, the virtual probe is in free motion, i.e., it is not in contact with the virtual surface.

• If ${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)=\mathbf{0}$ the virtual probe is in contact with the virtual surface and it stays over the surface only, staying in constrained motion.

• If ${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)<\mathbf{0}$ virtual probe is in restricted movement, but we violate restriction i.e., it is inside of the virtual surface.

For the virtual reality system, it is important to remember that the vector ${\mathit{x}}_{\mathrm{v}}$ represents the virtual robot’s end-effector position. The virtual probe would share such a position. By extending the approach proposed, the virtual robot acts in fact as the virtual tool and it projected its position to the operator by the avatar of the system. In the second case, for an nonpenetrable constraint, even when the contact starts at a single point, the properties of the surface allow the virtual probe to stay in motion, as we can see in Figure 7. Intuitively, the process is more complex since, at a certain moment, the virtual probe must stop. In a medical context, that means the deformable tissue has a limited resistance that depends on its elastic properties. Since the aim is to get a reaction force that depends on the motion of the virtual probe once inside the deformable object, it is necessary to describe such motion properly. Unlike the non-deformable surfaces, the force does not arise normally to the surface at a single point, but lateral forces occur when the operator tries to move the virtual probe in such directions. Ideally, this would be true for any method to represent soft tissues, including finite element meshes and particle-based models. Figure 8 shows the virtual tool motion inside a virtual object.

Figure 8.

Motion in 2D of the virtual tool inside a virtual object.

For easy visualization, it showed the motion in 2D but during the simulation, we must reproduce it in 3D with the aim of increase the realism of the application by improving the operator’s dexterity. The contact begins in stage A where the virtual tool penetrates the object by following a straight trajectory, represented by a blue arrow, to reach the position in stage B. It can follow other trajectories, represented by dashed red lines, to reach the position of stage D or stage C. However, because of the surrounding tissue, the tool cannot move laterally because of the reaction forces (represented by the red arrows) preventing it along the trajectories. In contrast to what happens in the holonomic case, the virtual tool may stay in motion, i.e., its velocity differs from zero until the operator stops voluntarily. The process described above is like the motion of a wheeled car in 2D, which is perfectly described by non-holonomic constraints. In fact, if we add a third dimension, Property 2.2 remains forever and the virtual tool can reach any point on the virtual object.

The trajectories of the virtual tool shown in Figure 8 are common in noninvasive surgical procedures. For example, in the simulator of Transurethral Resection of the Prostate (TURP Mentor) of Figure 2, the medical trainee performs a straight trajectory that, in the first place, simulates the insertion of the resectoscope into the patient’s penis. Once within the virtual prostate, the resectoscope needs to be moved to remove, through an incandescent resection loop, the benign tissue that obstructs the flow of urine to the urethra. These movements follow a similar route to the one represented with red dotted lines in Figure 8. These movements use a pivot where the tool changes its direction.

The contribution of this approach is, in contrast to common single point haptic feedback methods, that we produce forces that prevent lateral movement of the virtual tool. However, the major disadvantage is that the environment’s elastic properties are not considered. As a result, the reaction of the force that limits the movement of the operator, depending on such properties, is not reproduced and can move the virtual tool interchangeably within the virtual object, which does not happen in real life. For example, in human organs, elastic properties and parameters such as Young’s modulus or Poisson’s ratio establish limits of movement for the tool that, when exceeded, the tissue is damaged.

3.1.1 Virtual environment design

In Figure 9, we show a scheme in 2D of the contact between a virtual tool and a virtual surface to illustrate the use of the model given by Eq. (21).

Figure 9.

Virtual robot in interaction with a penetrable surface.

We attach the tool to the virtual robot’s last DOF, acting as its end-effector. It is important to note that the manipulator dynamic model is used to reproduce a classic bilateral teleoperation system and assuming that, since it is simulated digitally, it can be exchanged by a simpler or more complex model, including those of medical instruments such as forceps endoscopes, gripers, and retractors. This assumption leads to the proposition that, if we use the model of a surgical tool during the simulation, the realism of the contact with the surface would increase. The principal difference between defining a holonomic and a non-holonomic constraint is a need for an expression of ${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)$. Based on Section 2.1.1, from an implicit representation approach, we build rigid virtual objects from 3D basic geometric primitives as cones, pyramids, planes, cubes, and spheres, [22]. Ultimately, the base of a highly complex virtual environment composed of rigid objects is a set of basic geometric shapes that we can represent through mathematical expressions. Therefore, it is enough to define a zero set of functions as in (25) that are individually expressed as ${\mathit{\phi }}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)$ and a collision detection algorithm based on the inequalities stated above. For non-holonomic constraints and considering again the virtual robot of Figure 9, let ${}^0{\mathit{p}}_{\mathrm{v}}\in {\mathbb{R}}^3$ be the Cartesian position of the virtual robot end–effector and ${}^0{\mathit{R}}_{\mathrm{v}}\in {\mathit{SO}}^3$ a rotation matrix that describes its orientation. Dividing this rotation matrix into three column vectors as

$${}^0{\mathit{R}}_{\mathrm{v}}=\left[\begin{array}{ccc}{}^0{\mathit{x}}_{n\mathrm{v}}& {}^0{\mathit{y}}_{n\mathrm{v}}& {}^0{\mathit{z}}_{n\mathrm{v}},\end{array}\right]$$E38

for which each column represents a vector of the end-effector coordinate frame, described in the base frame. This allows defining Pfaffian constraints like (29) in an intuitive form, i.e.,

$${\mathit{A}}_{\mathrm{v}}\left({}^0{\mathit{x}}_{n\mathrm{v}}{,}^0{\mathit{y}}_{n\mathrm{v}}{,}^0{\mathit{z}}_{n\mathrm{v}}\right){\mathit{v}}_{\mathrm{v}}=\mathbf{0}.$$E39

We claim that a set of non-holonomic constraints can be defined if the manipulator degrees of freedom are greater than those necessary to control the end-effector position, i.e., $n>2$ for planar robots and $n>3$ for robots in a three-dimensional workspace. The end-effector velocities of the virtual robot can be described by

$${\mathit{v}}_{\mathrm{v}}=\left[\begin{array}{c}{}^0{\dot{\mathit{p}}}_{\mathrm{v}}\\ \\ {}^0{\mathit{\omega }}_{\mathrm{n}}\end{array}\right]$$E40

where ${\mathit{\omega }}_{\mathrm{n}}$ is the angular velocity over an axis normal to the robot plane, and ${}^0{\dot{\mathit{p}}}_{\mathrm{v}}$ is the linear velocity defined as

$${\dot{\mathit{p}}}_{\mathrm{v}}=\left[\begin{array}{c}{}^0{\dot{\mathit{p}}}_{\mathrm{v}x}\\ \\ {}^0{\mathit{p}}_{\mathrm{v}y}\end{array}\right]$$E41

If the robot may not move in the ${}^0{\mathit{y}}_{n\mathrm{v}}$ direction, the corresponding Pfaffian constraint is given by

$$\left[\begin{array}{cc}{}^0{\mathit{y}}_{n\mathrm{v}}^{\mathrm{T}}& 0\end{array}\right]{\mathit{v}}_{\mathrm{v}}=\left[\begin{array}{ccc}-sin\left(q_{\mathrm{v}1}+q_{\mathrm{v}2}+q_{\mathrm{v}3}\right)& cos\left(q_{\mathrm{v}1}+q_{\mathrm{v}2}+q_{\mathrm{v}3}\right)& 0\end{array}\right]\left[\begin{array}{c}{}^0{\dot{p}}_{\mathrm{v}x}\\ {}^0{\dot{p}}_{\mathrm{v}y}\\ {}^0{\dot{\omega }}_n\end{array}\right].$$E42

By choosing the distribution

$$\mathrm{\Delta }=\left[\begin{array}{c}{\mathit{g}}_1,{\mathit{g}}_2\end{array}\right]$$E43

where

$${\mathit{g}}_1=\left[\begin{array}{c}cos\left(q_{\mathrm{v}1}+q_{\mathrm{v}2}+q_{\mathrm{v}3}\right)\\ sin\left(q_{\mathrm{v}1}+q_{\mathrm{v}2}+q_{\mathrm{v}3}\right)\\ 0\end{array}\right]$$E44

$${\mathit{g}}_2=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$$E45

as a basis for the null-space of the Pfaffian matrix, the equivalent control system

$${\dot{\mathit{q}}}_{\mathrm{v}}={\mathit{g}}_1{u}_1+{\mathit{g}}_2{u}_2$$E46

can be constructed, representing the directions of allowed motion, [42].

It is easy to verify that the Lie bracket is

$$\left[\begin{array}{cc}{\mathit{g}}_1& {\mathit{g}}_2\end{array}\right]=\left[\begin{array}{c}-sin\left(q_{\mathrm{v}1}+q_{\mathrm{v}2}+q_{\mathrm{v}3}\right)\\ cos\left(q_{\mathrm{v}1}+q_{\mathrm{v}2}+q_{\mathrm{v}3}\right)\\ 0\end{array}\right],$$E47

which shows the non-involutivity of the distribution and thus establishes the non-holonomic nature of the constraints according to Proposition 2.1.

Notice also that if the degrees of freedom were 2, the null-space would be of dimension 1, which is necessarily involutive, and the constraints would be holonomic.

3.2 Position-force controllers design

A correct haptic rendering largely depends on the force control algorithm. In classic haptic systems, the common solution is to define indirect impedance or compliance control schemes. In contrast, in this section, we present two hybrid control algorithms for haptic interaction with virtual constrained systems. As mentioned in Section 2, the usual practice is to associate the position of the haptic interface directly with that of the virtual avatar. Therefore, a position control scheme is unnecessary, as the operator’s movements are reflected in the graphical application accurately. However, in the proposed approach, the task space coordinates of the virtual environment depend on the correct tracking between the position of the haptic robot and that of the virtual one, i.e., the control algorithm generates the virtual environment itself. This is due to including the virtual robot dynamics and the fact that the operator should feel the virtual tool because of the masking effect. In order to address this, we explored a control scheme used in teleoperation to achieve both position and force tracking. Next, we show a block diagram of this scheme in Figure 10.

Figure 10.

Block diagram of the proposed scheme.

Considering once again $i,j=\mathrm{m},\mathrm{v}$ where $i\ne j$. We define

$${\mathit{q}}_{\mathrm{d}i}\left(t\right)\triangleq {\mathit{q}}_j\left(t\right)$$E48

as the desired position trajectories, and

$${\dot{\mathit{q}}}_{\mathrm{d}i}\left(t\right)\triangleq {\dot{\mathit{q}}}_j\left(t\right)$$E49

as the desired velocity trajectories, i.e., if $i=\mathrm{m}$, then $j=\mathrm{v}$ and vice versa.

We define the corresponding tracking error as

$$\mathrm{\Delta }{\mathit{q}}_i\triangleq {\mathit{q}}_i-{\mathit{q}}_{\mathrm{d}i}.$$E50

Based on [48], we proposed

$${\mathit{s}}_i={\dot{\mathit{q}}}_i-{\dot{\mathit{q}}}_{\mathrm{d}i}+{\mathbf{\Lambda }}_{\mathit{xi}}\mathrm{\Delta }{\mathit{q}}_i.$$E51

and

$${\dot{\mathit{\sigma }}}_i={\mathit{K}}_{\mathit{\beta i}}{\mathit{s}}_i+sign\left({\mathit{s}}_i\right),$$E52

where ${\mathit{K}}_{\mathit{\beta i}}\in {\mathbb{R}}^{n\times n}$ is a positive definite diagonal matrix and

$$sign\left({\mathit{s}}_i\right)=\left[\begin{array}{c}sign\left({\mathit{s}}_{i1}\right)\\ ⋮\\ sign\left({\mathit{s}}_{\mathit{\text{in}}}\right)\end{array}\right]$$E53

with $s_{\mathit{ij}}$ element of ${\mathit{s}}_i$ for $j=1,\dots ,n$.

Now, by considering the velocity reference as

$${\dot{\mathit{q}}}_{\mathrm{r}i}={\dot{\mathit{q}}}_{\mathrm{r}i}+{\mathbf{\Lambda }}_{\mathit{xi}}\mathrm{\Delta }{\mathit{q}}_i-{\mathit{K}}_{\mathit{\gamma i}}{\mathit{\sigma }}_i,$$E54

where ${\mathit{K}}_{\mathit{\gamma i}}{\mathit{\sigma }}_i\in {\mathbb{R}}^{n\times n}$ is a positive definite diagonal matrix.

We define also the auxiliary variable

$${\mathit{s}}_{\mathrm{a}i}={\dot{\mathit{q}}}_i-{\dot{\mathit{q}}}_{\mathrm{r}i}.$$E55

By supposing that both robots are in free movement, for that case, the control laws for the master and the virtual robots are proposed as

$${\mathit{\tau }}_{\mathrm{m}}=-{\mathit{K}}_{\mathrm{am}}{\dot{\mathit{q}}}_{\mathrm{m}}-{\mathit{K}}_{\mathrm{pm}}{\mathit{s}}_{\mathrm{am}}$$E56

$${\mathit{\tau }}_{\mathrm{v}}={\mathit{K}}_{\mathrm{av}}{\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{K}}_{\mathrm{pv}}{\mathit{s}}_{\mathrm{av}},$$E57

respectively, where ${\mathit{K}}_{\mathrm{am}},{\mathit{K}}_{\mathrm{pm}},{\mathit{K}}_{\mathrm{av}}$ and ${\mathit{K}}_{\mathrm{pv}}$ are positive definite diagonal matrices.

3.2.1 Virtual holonomic constraints

Making an approximation of what happens during the tactile interaction of a point probe with a rigid surface, we considered the one-dimensional case (${\mathit{\phi }}_{\mathrm{v}}:{\mathbb{R}}^n\to \mathbb{R}$). As mentioned in Section 3.1, that is ideally the normal force generated at a single point of contact where other reactions, as friction or tangential forces, we can omit them, [32]. In order to reproduce this effect, we use the implicit surface method, which ${\lambda }_{\mathrm{v}}={\lambda }_{\mathrm{v}}$ represents the normal force of the virtual manipulator over the virtual surface. To reflect such contact force, the Generic Penalty Method computed a Lagrange multiplier as used by [49], i.e.,

$${\lambda }_{\mathrm{v}}={\alpha }_{\mathrm{v}}\left({\ddot{\mathit{\phi }}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)+2{\mathit{\xi \omega }}_{\mathrm{n}}{\dot{\mathit{\phi }}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)+{\omega }_n^2{\mathit{\phi }}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)\right)$$E58

where $\xi ,{\omega }_{\mathrm{n}}>0$. Considering that force measurements are available at the master side in Cartesian coordinates and mapping the virtual force to this space as

$${\mathit{F}}_{\mathrm{v}}={\mathit{J}}_{\mathit{\phi x}\mathrm{v}}^{\mathrm{T}}\left({\mathit{x}}_{\mathrm{v}}\right){\lambda }_{\mathrm{v}},$$E59

we use a PID-like controller for the virtual reality system. Consider that ${\mathit{F}}_{\mathrm{h}}\in {\mathbb{R}}^3$ is the normal force component measured with a force sensor mounted at the master robot end-effector. After (23), (24), and (58), we can use a PI controller for the virtual reality system.

By defining

$${\mathit{F}}_{\mathrm{d}i}\left(t\right)\triangleq {\mathit{F}}_i\left(t\right)$$E60

as the desired force trajectory where if $i=\mathrm{h}$ then $j=\mathrm{v}$ and vice versa, as stated before.

The force tracking errors are

$$\mathrm{\Delta }{\mathit{F}}_i={\mathit{F}}_i-{\mathit{F}}_{\mathrm{d}i}$$E61

and the corresponding integral, the momenta tracking error, is

$$\mathrm{\Delta }{\mathit{\rho }}_i={\int }_0^t\mathrm{\Delta }{\mathit{F}}_i\mathrm{dt}.$$E62

Note that we use the standard notation for momenta $\mathit{\rho }$, although also the same notation is for the position. We claim that there is no confusion because it always appears $\mathrm{\Delta }\mathit{\rho }$ in that case. Instead of using Eqs. (56) and (57) to describe the master and the virtual robot, we gave the corresponding control laws by

$${\mathit{\tau }}_{\mathrm{m}}={\mathit{Y}}_{\mathrm{m}}\left({\mathit{q}}_{\mathrm{m}},{\dot{\mathit{q}}}_{\mathrm{m}},{\ddot{\mathit{q}}}_{\mathrm{m}}\right){\mathit{\theta }}_{\mathrm{m}}-{\mathit{K}}_{\mathrm{am}}{\dot{\mathit{q}}}_{\mathrm{m}}-{\mathit{K}}_{\mathrm{pm}}{\mathit{s}}_{\mathrm{am}}+{\mathit{J}}_{\mathrm{m}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{m}}\right)\left({\mathit{F}}_{\mathrm{v}}-{\mathit{K}}_{\mathrm{fm}}\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{h}}\right)$$E63

$${\mathit{\tau }}_{\mathrm{v}}={\mathit{K}}_{\mathrm{av}}{\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{K}}_{\mathrm{pv}}{\mathit{s}}_{\mathrm{av}}-{\mathit{J}}_{\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right)\left({\mathit{F}}_{\mathrm{h}}-{\mathit{K}}_{\mathrm{fv}}\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{v}}\right),$$E64

respectively, where ${\mathit{K}}_{\mathrm{f}i}\in {\mathbb{R}}^{n\times n}$ are diagonal matrices. In order for the operator to feel the virtual tool in contact with the virtual environment, we can carry out a dynamic cancelation of the master manipulator dynamics, as shown in (63).

3.2.2 Virtual non-holonomic constraints

In contrast with the holonomic case, when the constraints are non-holonomic, we cannot define them as a function of a set of generalized coordinates, as stated by the Frobenius theorem. As a result, we cannot compute the Lagrange multipliers as in (58). We define these constraints in the form (31) or equivalently (34). One problem arising from these constraints is how to compute the Lagrangian multipliers to satisfy (36). These multipliers represent the forces required to maintain such constraints. Unfortunately, most of the methods used to calculate the lagrange multipliers are designed for systems with holonomic constraints [30, 49, 50] and, therefore, these methods require a position-level definition of the Pfaffian constraints as in (25) or (27). As stated in (45), the calculation presented in [42] can be used for this case. However, it is well known that this solution is unstable since its underlying mechanism is a second-order integrator with zero input. In this work, a modification of the approach used in [50] is proposed as follows. For simplicity, we define

$${\mathit{H}}_{\mathrm{v}}={\mathit{H}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)$$E65

$${\mathit{C}}_{\mathrm{v}}={\mathit{C}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}},{\dot{\mathit{q}}}_{\mathrm{v}}\right)$$E66

$${\mathit{g}}_{\mathrm{v}}={\mathit{g}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)$$E67

$${\mathit{A}}_{\mathrm{v}}={\mathit{A}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right)$$E68

$${\mathit{\psi }}_{\mathrm{v}}={\mathit{\psi }}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}},{\dot{\mathit{q}}}_{\mathrm{v}}\right)={\mathit{A}}_{\mathrm{v}}\left({\mathit{q}}_{\mathrm{v}}\right){\dot{\mathit{q}}}_{\mathrm{v}}$$E69

Then, the Lagrange multipliers can be computed as

$${\mathit{\lambda }}_{\mathrm{v}}={\left({\mathit{A}}_{\mathrm{v}}{\mathit{H}}_{\mathrm{v}}^{-1}{\mathit{A}}_{\mathrm{v}}^{\mathrm{T}}\right)}^{-1}\left[\dot{\mathit{\psi }}-{\dot{\mathit{A}}}_{\mathrm{v}}{\dot{\mathit{q}}}_{\mathrm{v}}-{\mathit{A}}_{\mathrm{v}}{\mathit{H}}_{\mathrm{v}}^{-1}\left({\mathit{\tau }}_{\mathrm{v}}-{\mathit{C}}_{\mathrm{v}}{\dot{\mathit{q}}}_{\mathrm{v}}-{\mathit{D}}_{\mathrm{v}}{\dot{\mathit{q}}}_{\mathrm{v}}-{\mathit{g}}_{\mathrm{v}}\right)\right],$$E70

where the constraints are forced to satisfy

$$\dot{\mathit{\psi }}+2{\alpha }_{\mathrm{v}}{\mathit{\psi }}_{\mathrm{v}}+{\beta }_{\mathrm{v}}{\int }_{t_0}^t{\mathit{\psi }}_{\mathrm{v}}\mathrm{d}\vartheta =\mathbf{0}$$E71

with ${\alpha }_{\mathrm{v}},{\beta }_{\mathrm{v}}>0$ chosen to ensure rapid convergence to the origin.

Note that the constraint function $\mathit{\psi }$ can be defined in terms of the velocities of the end-effector,

$$\mathit{\psi }=\mathit{\psi }\left({\mathit{x}}_{\mathrm{v}},{\mathit{v}}_{\mathrm{v}}\right)={\mathit{A}}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right){\mathit{v}}_{\mathrm{v}}.$$E72

Therefore, the initial condition of the integral term on the left-hand side of (72) can be set to zero. Each element ${\mathit{\lambda }}_{\mathrm{v}}$ is a function of ${\mathit{q}}_{\mathrm{v}},{\dot{\mathit{q}}}_{\mathrm{v}}$ and ${\mathit{\tau }}_{\mathrm{v}}$ since the constraints change with the configuration, velocity and virtual applied force.

By substituting (71) in the motion Eq. (72), a complete description of the dynamics of the system is gotten. Regarding force, sensor measurements ${\mathit{F}}_{\mathrm{v}}$ on the master side can calculate the real Lagrange multiplier as

$${\mathit{\lambda }}_{\mathrm{m}}={\left({\mathit{A}}_{\mathrm{v}}{\mathit{A}}_{\mathrm{v}}^{\mathrm{T}}\right)}^{-1}{\mathit{A}}_{\mathrm{v}}{\mathit{F}}_{\mathrm{h}}.$$E73

By defining

$${\mathit{\lambda }}_{\mathrm{d}i}\left(t\right)\triangleq {\mathit{\lambda }}_{\mathrm{j}}\left(t\right)$$E74

as the desired force trajectory in joint space. The corresponding integral is

$$\mathrm{\Delta }{\mathit{\lambda }}_i={\int }_0^t\left({\mathit{\lambda }}_i-{\mathit{\lambda }}_{\mathrm{d}j}\mathrm{dt}\right)$$E75

Finally, instead of (56) and (57), for the master and virtual robot the proposed position-force control for a virtual dynamic system subject to non-holonomic constraints is

$${\mathit{\tau }}_{\mathrm{m}}={\mathit{Y}}_{\mathrm{m}}\left({\mathit{q}}_{\mathrm{m}},{\dot{\mathit{q}}}_{\mathrm{m}},{\ddot{\mathit{q}}}_{\mathrm{m}}\right){\mathit{\theta }}_{\mathrm{m}}-{\mathit{K}}_{\mathrm{am}}{\dot{\mathit{q}}}_{\mathrm{m}}-{\mathit{K}}_{\mathrm{pm}}{\mathit{s}}_{\mathrm{am}}+{\mathit{A}}_{\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right)\left({\mathit{\lambda }}_{\mathrm{v}}-{\mathit{K}}_{\mathrm{fm}}\mathrm{\Delta }{\mathit{\lambda }}_{\mathrm{m}}\right)$$E76

$${\mathit{\tau }}_{\mathrm{vm}}={\mathit{K}}_{\mathrm{av}}{\dot{\mathit{q}}}_{\mathrm{v}}+{\mathit{K}}_{\mathrm{pv}}{\mathit{s}}_{\mathrm{av}}-{\mathit{A}}_{\mathrm{v}}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{v}}\right)\left({\mathit{\lambda }}_{\mathrm{m}}-{\mathit{K}}_{\mathrm{fv}}\mathrm{\Delta }{\mathit{\lambda }}_{\mathrm{v}}\right)$$E77

Note that the novelty of the approach is not the control scheme because very well-known techniques are employed, but the novelty lies in the effective use of non-holonomic constraints to describe penetrable virtual surfaces. Therefore, a technical stability proof is not provided, but it shows a set of reliable experiments in the next section with the aim of validating the proposed approach.

3.3 Visual components of the virtual environment

A fundamental part of the developed virtual reality system is visual feedback. In dynamic systems and control research, there is no interest in including such elements but in real-world applications, as surgery simulators, it is essential. Nowadays, in those developments, we compose the virtual environments by merging several numeric techniques that, combined with the fast velocity of today’s processors, give the virtual objects and surfaces a realism that before would seem impossible to reach., since the goal of this dissertation is to show how a teleoperation control scheme can be used in a virtual reality system, we design the environment by using the fundamentals of graphic computing. The tool used to design the virtual environment was the graphic standard OpenGL 2.0 which is an API, which is a software library for accessing features in graphics hardware. It contains different commands that are used to specify objects, images, and operations needed to produce interactive three-dimensional graphic applications, [51]. Among those operations, the possibility to give texture¹ and lighting to the virtual objects is possible, besides proportioning position and orientation changes to the scene’s camera, i.e., the way the operator sees the images on the computer screen regarding height, deep, viewing angles as pitch, roll, and yaw, etc. As we can see in Figure 11, the environment of the developed application comprises a motionless floating sphere and the virtual avatar of the system.

Figure 11.

Virtual environment developed in OpenGL 2.0.

For simplicity’s sake, there are no changes in the camera’s position and orientation, but we gave lighting and texture to the scene. We appreciate a notable difference in Figure 12, where the avatar has no lighting, the quality of its texture is less and the background color changes.

Figure 12.

Virtual avatar of the system.

Here, a set of aligned cylinders becomes the avatar, which directly related position and orientation to those of the end-effector of the haptic interface. Since OpenGL has the instructions to create elements from primitives, the generation of such cylinders and the floating sphere was straightforward.

3.3.1 Rigid sphere

It is important to remember that, from a haptic rendering approach, we established the classification of nonpenetrable and penetrable objects for rigid and deformable objects, respectively. With the rigid object, the Open GL instruction glutSolidSphere() build automatically a solid sphere with a specific radius by defining the number of subdivisions around (sectors) and along (stacks) the $z$ axis, as we can see in Figure 13.

Figure 13.

Sectors and stacks of a solid sphere.

We give the effect of rigidity because that vertex’s position is not changed when contact with the virtual avatar occurs. However, giving the haptic effect of highly rigid objects to the operator was difficult since an impedance device was used. For such reason, it is important to establish the control scheme (63) and (64) that compensates, as possible, the limitations of hardware and make feel it produced contact with a rigid object.

3.3.2 Deformable sphere

The real challenge comes when the virtual object is deformable. Saying that an object is deformable has many implications, mainly related to mechanics. Therefore, the visual effect of deformation is more complex than that produced by stiffness, because a real deformable object has an infinite number of degrees of freedom. For this reason, virtual objects need a high resolution, which gives a better rendering quality to visual and haptic feedback [52] and more realism to the application. However, this is always limited to the computational resources available. We must run both the graphics and the control algorithm in a single program, with the smallest sample time that the system allows. If the graphics part occupies more processing resources, this sampling time will increase and there will be unwanted effects, delays, and, finally, an application crash. We drew the sphere by defining four vertices $a,b,c$ and $d$, which form a plane that is replicated iteratively according to several parallels $p$ and meridians $m$ defined by the operator. We intrinsically linked the value of the iteration to resolving the sphere. Figure 14, it showed the deformable sphere with different resolutions. In the sphere on the left, the values used were $p=m=50$ while in the central sphere were $p=m=100$. For the sphere on the right, the values were $p=m=150$.

Figure 14.

Surface mesh generated using different values for $p$ and $m$. (a) p = m = 50. (b) p = m = 100. (c) p = m = 150.

As mentioned before, the more the resolution of the sphere, the more the realism of the application. However, the computational processing when using the resolution of the last case did not allow a correct performance of the graphic or haptic part. For this reason, a resolution $p=m=100$ was chosen. For simplicity’s sake, the sphere was built by placing two hemispheres, one above the other, regarding a common axis. The contact with the virtual avatar will arise in a single point $\left(x_{\mathrm{v}},y_{\mathrm{v}},z_{\mathrm{v}}\right)$ computed parametrically as

$$x_{\mathrm{v}}=rcos\left(\alpha \right)cos\left(\beta \right)$$E78

$$y_{\mathrm{v}}=rcos\left(\alpha \right)sin\left(\beta \right)$$E79

$$z_{\mathrm{v}}=rsin\left(\alpha \right)$$E80

where $r$ is the radius of the sphere, and $\alpha$ and $\beta$ are the angles from whose ranges parallel and meridians are drawn. For code optimization, meridians have the range of $\left[\mathrm{0,360}\right]$ and parallels of $\left[\mathrm{0,180}\right]$. Such ranges correspond to the upper hemisphere while the lower is drawn by considering the $z_{\mathrm{v}}$ axis negative part.

Every vertex $a,b,c$ and $d$ must take the value of Eqs. (78)–(80) in order to visualize its initial position in the virtual environment. To improve the interaction with the virtual avatar, we include an offset $r_{\mathrm{off}}$ in the equations that define the vertices as

$$v_{\mathit{nx}}=\left(r-r_{\mathrm{off}}\right)cos\left(\alpha \right)cos\left(\beta \right)$$E81

$$v_{\mathit{ny}}=\left(r-r_{\mathrm{off}}\right)cos\left(\alpha \right)sin\left(\beta \right)$$E82

$$v_{\mathit{nz}}=\left(r-r_{\mathrm{off}}\right)sin\left(\alpha \right)$$E83

where $n=a,b,c,d$ and $r_{\mathrm{off}}$ take an arbitrary value defined experimentally. The contact with the virtual avatar occurs in some plane defined iteratively by (81)–(83) using a collision detection algorithm consisting on validating the value of each vertex of every plane and comparing the value of each component $\left(x_{\mathrm{v}},y_{\mathrm{v}},z_{\mathrm{v}}\right)$ with the position of the master robot ${\mathit{p}}_{\mathrm{m}}$. If such values belong to the range of the plane, then the avatar is in contact with the sphere. The next step is to produce the effect of motion of the contact plane. Algorithm 1 shows the pseudocode of such a process, which uses an auxiliary normal force in which direction the plane will move. This force does not belong to that calculated using the non-holonomic constraint but the gradient of Eq. (29) and we use it only for visual effect. The process described above makes only one plane to move and, if that happens, the deformation effect is not realistic. For this reason, moving the adjacent planes is necessary. It moved the more planes, the more realistic the effect of the object deformation. However, in the application developed, we changed only the position of the surrounding planes to the contact plane, since the method is purely geometric and not based on continuum mechanics. Using such an approach implies a lot of considerations that are beyond this research.

Algorithm 1 shows the part of the pseudocode where the modification of the adjacent planes takes place. This does not occur at the same time as the modification of the contact plane, since the code has not yet created these. Figure 15 shows the visual effect that the motion, both of the contact plane and adjacent planes produce. It is important to note that the code implemented needs to be optimized and, above all, to be adapted to the force rendering algorithm through non-holonomic constraints.

Figure 15.

Deformation effect of the sphere. (a) Contact. (b) Deformation effect.

4. Experimental platform

The experimental platform comprises a Geomagic Touch haptic robot with six revolute joints, having only the first three of them actuated. An ATI Nano–17 six-axis force sensor is adapted at the last link, as shown in Figure 16. A PC executes the control loop with a sample time of $T=2$ milliseconds.

Figure 16.

Experimental platform.

As mentioned in Section 3, the virtual environment comprises a sphere developed using the graphic standard OpenGL 2.0. We should note that both the control algorithm and the graphic simulation run in the same application developed in Visual Studio/C++.

A practical limitation of the Geomagic Touch robot is that it actuated only the first three joints. Therefore, a projection of both the force reflection and the controller torques is necessary, i.e., the contribution of the last two joint torques is neglected. The virtual robot does not have this limitation, and therefore is considered to be fully actuated. The master robot limitation is not so restrictive, since the virtual environment considers only force but not end-effector torque feedback, avoiding the problem of sensor/actuators asymmetry haptic interfaces [53]. The contribution of the last two joints to the force reflection is much less in magnitude when compared with the contribution of the first three joints.

4.1 Task description

A detailed description of the interaction process between the virtual tool and the virtual environment is presented, simulating separately a rigid and a penetrable sphere. Since the goal of this research was to extend the use of the control scheme to medical training applications, we adopted the shape of the avatar as a needle, as we can see in Figure 11. In medicine, procedures that use this tool are very common, with needle insertion being the most studied and simulated procedure [54]. In this procedure, the operator takes a sterile needle and slowly brings it closer to the patient, once in contact, the operator must be very careful and, through a tactile sensation, know if the soft tissue (muscle, organ) or rigid (bone) has been affected. In both cases, the contact surface produces a reaction force in opposition to the operator’s movements.

While on a rigid surface, the force does not let the needle penetrate the tissue, in a penetrable surface this is possible. The force behaviors are different, as in the first case, there is a major contribution in the normal direction, which would allow the operator to move the needle laterally over the surface. In the second case, the normal force contribution is smaller and the surrounding tissue would not allow moving the needle in the lateral directions.

In the approach presented, we assume we attach the virtual tool to the end-effector of a five degrees-of-freedom manipulator, which is not visible in the graphic simulation. It may seem counterintuitive because, evidently in real life, a needle does not have such dynamics. We use the robot model as a demonstrative example of other medical tools such as an endoscope, resectoscope, forceps. Attached to teleoperated surgical robot arm have such complex dynamics that must be modeled. The graphic simulations in those cases include pulling, gripping, clamping, and cutting, and therefore it is convenient to have a complete description of both the kinematics and the dynamics of the tool-tissue interaction, [3]. The task starts with the Geomagic Touch robot in its home position. The operator grasps the master robot stylus using the force sensor adapter tip to later gently bring it closer to the virtual surface. We imposed the desired trajectory in free motion in this way. The virtual robot moves following such a trajectory in the virtual environment, with no scaling between the virtual and the real workspaces. It perceived visually both the avatar movement and the virtual surface through a computer screen. When the collision-detection algorithm detects contact with the surface, the force-response algorithm generates a virtual force trajectory by computing the corresponding Lagrangian multipliers, either by employing (58) or (70). The operator perceives an interaction force exerted by the master robot and registered by the Nano-17 force sensor until the contact is over. Finally, the operator returns the sensor adapter to its initial position, thus completing the task.

4.2 Holonomic constraint experiment

For simplicity’s sake, the surface used to test the validity of the proposed approach is a sphere described by

$${\phi }_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right)={\left(x_{\mathrm{v}}-h\right)}^2+{\left(y_{\mathrm{v}}-k\right)}^2+{\left(z_{\mathrm{v}}-l\right)}^2-r^2=0,$$E84

where ${\mathit{x}}_{\mathrm{v}}={\left[\begin{array}{ccc}{x}_{\mathrm{v}}& y_{\mathrm{v}}& z_{\mathrm{v}}\end{array}\right]}^{\mathrm{T}}$ is the vector that stands for the virtual environment task-space coordinates, $r=0.1$ [m] is the radius, and $\left(h,k,l\right)=\left(\mathrm{0.4,0,0}\right)$ [m] are the sphere center coordinates. It is important to note that, in contrast with other works, we added a third dimension ${\mathbf{z}}_{\mathrm{v}}$ in order to heighten the realism of the virtual reality application. For example, in [44, 49], the authors consider only two dimensions to test different control schemes for a haptic and a teleoperation system, respectively. The gains of the master manipulator, described in (63), and the virtual manipulator, are shown in Table 1.

Variable (control law)	Value	Variable (virtual manipulator)	Value
${\mathit{K}}_{\mathrm{am}}$	0.0550 $\mathit{I}$	${\mathit{K}}_{\mathrm{av}}$	0.20 $\mathit{I}$
${\mathit{K}}_{\mathrm{pm}}$	0.0055 $\mathit{I}$	${\mathit{K}}_{\mathrm{pv}}$	$diag\left\{\mathrm{0.2,0.2,0.2,0.1,0.1}\right\}$
${\mathit{K}}_{\mathrm{fm}}$	10.050 $\mathit{I}$	${\mathit{K}}_{\mathrm{fv}}$	0.20 $\mathit{I}$
${\mathbf{\Lambda }}_{\mathrm{xm}}$	0.2500 $\mathit{I}$	${\mathbf{\Lambda }}_{\mathrm{xv}}$	20.0 $\mathit{I}$
${\mathit{K}}_{\beta \mathrm{m}}$	0.0100 $\mathit{I}$	${\mathit{K}}_{\beta \mathrm{v}}$	1.00 $\mathit{I}$
${\mathit{K}}_{\gamma \mathrm{m}}$	0.0150 $\mathit{I}$	${\mathit{K}}_{\gamma \mathrm{v}}$	0.20 $\mathit{I}$

Table 1.

Gains from the control law and the virtual manipulator (Holonomic constraint experiment).

$\mathit{I}$ is the identity matrix, which has the appropriate dimensions.

Finally, by using the Generic Penalty Method, the surface parameters are ${\alpha }_{\mathrm{v}}=0.002,\xi =100$ and ${\omega }_{\mathrm{n}}=200$.

4.3 Non-holonomic constraint experiment

As mentioned before, we cannot express a deformable surface implicitly, even when the operator perceives it as a sphere both visually and haptically. We use a discrete representation similar to that presented in [52], where we assume the surface to comprise many neighboring planes defined by shared nodes. We propose a technique that comprises iteratively choosing a small neighborhood of planes where the contact will occur, depending on the position of the virtual tool. Subsequently, we associate the Lagrange multipliers in Eq. (70) with a pair of planes using the impulse-based technique for multiple rigid body simulations, [55]. The micro-collisions with this technique occur only in the chosen vicinity of the sphere and Lagrange multipliers replace the impulses preventing body interpenetration. In the collision’s case detection algorithm, the implicit surface representation replaces the convex polyhedra decomposition, [56], using the Eq. (84). We did this for ease and to reduce the computational cost of the application, otherwise, the control algorithm sample time would increase. Considering the case where the needle is inside the sphere, but it may not move laterally. However, it may pivot to change orientation.

We adequately describe this kind of scenario by employing non-holonomic constraints. As mentioned in Section 1, non-holonomic constraints have been little exploited to represent the interaction with penetrable surfaces. For example, in [57] it is claimed that to model a surgeon’s scalpel both holonomic and non-holonomic constraints could be employed by limiting the depth of its incision and the direction of its motion, respectively. However, there is not any analysis or modeling of this process in such work.

The most representative proposal is the derivation of the non-holonomic generalized unicycle model presented in [58], where a coordinate-free representation is used to model the insertion of a flexible needle into soft tissue. We employed a similar approach by using the homogeneous matrix representation, but taking into consideration both the kinematic and the dynamic model of the virtual robot and the fact that non-holonomic constraints are more intuitively got if we define them in task space coordinates. The computation of the Lagrangian multipliers for non-holonomic constraints, which is proposed in (71), is an important improvement regarding the cited works. The experiment comprised five degrees of freedom virtual manipulator interacting with a deformable sphere. Once in contact, the end-effector may not move laterally, i.e., along the ${}^0{\mathit{y}}_{5\mathrm{v}}$ and ${}^0{\mathit{z}}_{5\mathrm{v}}$ axes, after a conventional Denavit-Hartenberg allocation, but it may move along the ${}^0{\mathit{x}}_{5\mathrm{v}}$ axis, i.e., along the pointing direction of the end-effector. The end-effector may rotate (pivoting) to change direction (as a three-dimensional version of the non-holonomic unicycle) and the Pfaffian matrix is computed as

$${\mathit{A}}_{\mathrm{v}}\left({\mathit{x}}_{\mathrm{v}}\right){\mathit{v}}_{\mathrm{v}}=\left[\begin{array}{cc}{}^0{\mathit{y}}_{5\mathrm{v}}& {\mathbf{0}}_{1\times 3}\\ {}^0{\mathit{z}}_{5\mathrm{v}}& {\mathbf{0}}_{1\times 3}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{\rho }}}_{\mathrm{v}}\\ {\mathit{\omega }}_{\mathrm{v}}\end{array}\right]=\mathbf{0}$$E85

The experiment has four steps, as shown in Figure 17. First, the virtual robot is in free motion and only the teleoperation part of the controller (77) is active. In the second part, we insert the needle into the sphere. Next, in the third part, approximately 45 degrees rotate the needle without changing its position. Finally, we insert the needle deeper into the sphere with a new orientation.

Figure 17.

Interaction sequence between the avatar and the non-holonomic virtual surface. (a) Free motion. (b) Insertion. (c) Pivoting. (d) Insertion.

The force of the human operator in the lateral directions of the needle is difficult to measure directly with the force sensor. Instead, we take advantage of the projection of such forces on the torque of the master manipulator joint, we calculate ${\mathit{\lambda }}_{\mathrm{m}}$ from (20), (23), and (73).

The gains of the master manipulator, described in (76), and the virtual manipulator, are shown in Table 2.

Variable (control law)	Value	Variable (virtual manipulator)	Value
${\mathit{K}}_{\mathrm{am}}$	0.0550 $\mathit{I}$	${\mathit{K}}_{\mathrm{av}}$	0.20 $\mathit{I}$
${\mathit{K}}_{\mathrm{pm}}$	0.0550 $\mathit{I}$	${\mathit{K}}_{\mathrm{pv}}$	$diag\left\{\mathrm{0.2,0.2,0.2,0.1,0.1}\right\}$
${\mathit{K}}_{\mathrm{fm}}$	0.010 $\mathit{I}$	${\mathit{K}}_{\mathrm{fv}}$	2.00 $\mathit{I}$
${\mathbf{\Lambda }}_{\mathrm{xm}}$	0.2500 $\mathit{I}$	${\mathbf{\Lambda }}_{\mathrm{xv}}$	20.0 $\mathit{I}$
${\mathit{K}}_{\beta \mathrm{m}}$	0.0150 $\mathit{I}$	${\mathit{K}}_{\beta \mathrm{v}}$	1.00 $\mathit{I}$
${\mathit{K}}_{\gamma \mathrm{m}}$	0.0150 $\mathit{I}$	${\mathit{K}}_{\gamma \mathrm{v}}$	0.20 $\mathit{I}$

Table 2.

Gains from the control law and the virtual manipulator (Non-holonomic constraint experiment).

$\mathit{I}$ is the identity matrix, which has the appropriate dimensions.

5. Conclusions

We present a proposal on chaptic interaction with holonomic and non-holonomic virtual constraints. Since extensive research on haptic interaction with rigid surfaces has been presented in the literature, the principal aim was to reproduce the forces generated by the interaction with soft surfaces from a force feedback approach.

Throughout the document, we introduced the theory to establish an optimal relationship between the visual and haptic interaction for the virtual reality system developed. The key lies in adapting the mathematical properties of the holonomic and non-holonomic constraints to the tool’s kinematics in contact with a nonpenetrable and penetrable virtual surface, respectively. However, it is important to note that we made this adaptation to achieve haptic feedback purposes only and consider the basic contact properties of simulated rigid and soft tissues.

Adapting a teleoperation control scheme to a virtual reality system was the strategy to follow since it allowed embedding a robot’s dynamic model into the virtual environment. By doing this, we addressed the teleoperated slave system as a problem of a virtual robot in constrained motion and either holonomic or non-holonomic constraints gave whose contact force.

We studied the differences between one or another representation, both mathematically and intuitively, and the particularities of each one. Among them, there is the fact that they could render forces using the Generic Penalty Method or the Pfaffian constraints matrix, respectively. We centered the interest on reproducing the contact with a soft surface by employing non-holonomic constraints.

The principal use detected is that such a method can render similar forces to those arising from contact between a tool and a penetrable surface. We can see, non-holonomic constraints have been used to reproduce the operator’s tactile sensations with practical meaning. We can eventually use this approach in virtual reality medical simulators, and it presented the fundamentals to do that throughout this document. However, adapting the developed method to complex virtual environments, such as those found in the medical field, requires more research both in control and computer graphics. In the first case, adapting more accurately the teleoperation controller presented is necessary and makes it fit with the current methods of graphic computing. Finding more optimal ways to model force using non-holonomic constraints is essential to heighten the realism of the applications. Regarding graphic computing, it is necessary to design numerical methods that adapt more efficiently to the control algorithms designed and that are capable of running with continuum mechanics models.

Naturally, this will increase computational processing and require more analysis to establish the compensation between real-time processing and control performance, without sacrificing application realism. All this requires a wide range of knowledge that does not belong to control or computer graphics.

Acknowledgments

The authors are grateful for the support provided by PRODEP-BUAP ID90934 and the DGAPA-UNAM under grant IN117820.

Algorithm 1: Deformable sphere

Require: Number of Parallels $p$.

Require: Number of Meridians $m$

Require: Radius $r$

Require: Offset $r_{\mathrm{off}}$

Define: Vertex $v_{\mathit{ax}},v_{\mathit{ay}},v_{\mathit{az}}$

Define: Vertex $v_{\mathit{bx}},v_{\mathit{by}},v_{\mathit{bz}}$

Define: Vertex $v_{\mathit{cx}},v_{\mathit{cy}},v_{\mathit{cz}}$

Evaluate:${\mathrm{\Delta }}_1={180}^{\circ }/p$ and ${\mathrm{\Delta }}_2={360}^{\circ }/m$

1: for all$i=0$ to $p/2$do

2: for all$j=0$ to $m$do

3: $\alpha =i\times {\mathrm{\Delta }}_1$

4: $\beta =j\times {\mathrm{\Delta }}_2$

5: Calculate the vertices of the upper hemisphere using (81)-(83).

6: if avatar touches a superior plane then

7: Move the contact plane through the auxiliary normal force

8: Store data of the contact plane (based on $p$ and $m$)

9: end if

10: Compute vertices of the lower hemisphere

11: if avatar touches an inferior plane then

12: Move the contact plane using (81)-(83) and the negative numbers on the $z_{\mathrm{v}}$ axis.

13: Store data of the contact plane (based on $p$ and $m$)

14: end if

15: end for

16: end for

17: if avatar touches a plane then

18: Reordering adjacent planes

19: end if

20: for all$i=0$ to $p/2$do

21: for all$j=0$ to $m$do

22: Draw plane

23: end for

24: end for