**Complete Geometric Analysis Using the Study SE(3) Parameters for a Novel, Minimally Invasive Robot Used in Liver Cancer Treatment**

**Iosif Birlescu <sup>1</sup> , Manfred Husty 2, Calin Vaida 1, Nicolae Plitea 1, Abhilash Nayak <sup>3</sup> and Doina Pisla 1,\***


Received: 8 November 2019; Accepted: 5 December 2019; Published: 7 December 2019

**Abstract:** The paper presents a complete geometric analysis of a novel parallel medical robotic system designed for minimally invasive treatment of hepatic tumors using brachytherapy, ablation or targeted chemotherapy. An algebraic method based on the study parameters of the special Euclidean transformation Lie group SE(3) was used to determine the mechanism kinematics singularities and workspace. Moreover, two particular medical tool manipulations for the minimally invasive medical procedures are defined in terms of the Study parameters. The first manipulation of the medical tool refers to the linear insertion (of e.g., needles) and the second one is the remote center of motion manipulation of specific medical instruments (e.g., ultrasound probes). The constraint equations of the robotic system are derived and then, the operational workspace is illustrated for the novel parallel robotic system. Lastly, a numerical simulation is presented showing the behavior of the robotic system manipulating the ultrasound probe constrained by the remote center of motion. The geometric analysis of the operational workspace and the numerical simulation show promising results that validate the novel robotic system (safe-wise) for the medical procedure.

**Keywords:** parallel robot; minimally invasive procedures; algebraic modeling; Study parameters

#### **1. Introduction**

In robotics and more specifically in mechanism analysis and synthesis there exists various mathematical approaches to describe the mechanism kinematics, singularities and workspace, which have both advantages and disadvantages [1]. Usually, the mathematical formulation used in describing a mechanism is chosen to benefit from the specific advantages that the method provides and depends strictly on the type (and even complexity) of the mechanism and sometimes on the task of the robot. One algebraic method based on the Study parameters (or dual quaternion) of SE(3) which is presented in [2–4], has the advantages (over some vector-based methods) that it describes the global kinematics of the mechanisms (showing all the working modes and being a powerful tool to find the mechanism singularities) and it is free of parameterization singularities. The Study parameters method was used in the analysis of various robotic systems such as the Stewart-Gough platform [5], the 3-RPS parallel manipulator [3], in the analysis of medical robots [4,6,7] and it was even tailored to analyze rehabilitation robots such as the one found in [8].

In this paper, the authors exploited the advantages of the algebraic method based on the Study parameters to achieve two main objectives: (i) a complete geometric analysis of the Pro-Hep-LCT robotic

system (which is designed for the liver cancer treatment through brachytherapy, ablation or targeted chemotherapy under ultrasound imaging), and (ii) the analysis of the behavior of the robotic system while it performs the medical task. The motivation of using the Study parameters for the first objective is simply due to the fact that the Study parameters offer a description of the global kinematics of the mechanism which in turn shows all the kinematic solutions, all the singularities and all the working modes of the robotic system. All this information, especially the singularities of the mechanism (since they define configurations where the control of the robotic system may be lost), should be valuable for understanding the robotic system capabilities to perform the medical task from the safety point of view. The motivation of the second objective was to provide the general Study parameter relations for two fundamental ways in which the medical instruments are manipulated in minimally invasive procedures (the linear insertion and the Remote Center of Motion (RCM) manipulation [9]) and use these general relations (in combination with the constraint equations of the Pro-Hep-LCT robotic system) to assess the robotic system behavior during the medical task (while the instruments are inserted inside the patient body) to determine the geometry of the operational workspace of the robot). The constraint equations (derived from the Study kinematic model) together with the general formulations for the medical insertion and RCM manipulation will enable the development of a robust fail-safe control algorithm. As explained in Ref. [10], interpolation may be used in combination to inverse kinematics of the robotic system to determine accurate relations between the medical tool coordinates (constrained by the RCM) and the active joints of the robot (this method was also demonstrated in Ref. [11] for serial robots). Moreover, in Ref. [12] the authors used a so called RCM-constrained Jacobian for the control of 7 DOFs and 9 DOFs MIS robotic systems. In Ref. [13] the authors described a collaborative framework for a 7 DOFs serial robot for minimally invasive surgery, showing also numerical simulations which describe the operating workspace of the robot. Line symmetric motion generators (which are sub manifold on the Study quadric) were presented in Ref. [14], showing their application in applications that require RCM.

The paper is structured as follows: Section 2 presents the mathematical framework of the Study parameters and it defines the mathematical relations for the linear insertion/retraction of a medical instrument and for the RCM manipulation; Section 3 presents the Pro-Hep-LCT robotic system and shows the kinematics, the singularity analysis and the workspace of the robotic system; Section 4 presents some numerical simulations for the Pro-Hep-LCT robotic system based on the previous obtained results; Section 5 presents the conclusions and future work.

#### **2. Mathematical Framework**

In the fields of multi-body systems kinematics (e.g., mechanism analysis), the special Euclidean group SE(3) describes every Euclidean displacement (with 6 DOFs) by means of homogeneous transformation matrices which encode both translations and orientation information. A kinematic mapping *k*, as found extensively in the scientific literature (see for example [2–4]), is defined as:

$$\begin{array}{c} k: SE(3) \to Q \in P^{\mathbb{Z}}\\ D(\mathbf{x}\_{i}, \mathbf{y}\_{i}) \to \begin{bmatrix} \mathbf{x}\_{0}: \mathbf{x}\_{1}: \mathbf{x}\_{2}: \mathbf{x}\_{3}: \mathbf{y}\_{0}: \mathbf{y}\_{1}: \mathbf{y}\_{2}: \mathbf{y}\_{3} \end{bmatrix}^{T} \neq \begin{bmatrix} 0: 0: 0: 0: 0: 0: 0: 0: 0 \end{bmatrix}^{T} \end{array} \tag{1}$$

and maps every Euclidean displacement *D* into a point *Q* from the projective space *P*7. The coordinates of the point *Q* are based on a dual quaternion and are also called the Study parameters in the scientific literature. The *xi* parameters show the orientation of a mobile coordinate frame (relative to a fixed one) whereas the parameters *yi* represent the translation between the origins of the two coordinate frames. The trivial solution (where all Study parameters are 0) is excluded since it has no use in kinematics. Moreover, the Study parameters must fulfill two conditions [2–4]:

$$\mathbf{x}\_0^2 + \mathbf{x}\_1^2 + \mathbf{x}\_2^2 + \mathbf{x}\_3^2 = 1 \tag{2}$$

$$x\_0 y\_0 + x\_1 y\_1 + x\_2 y\_2 + x\_3 y\_3 = 0\tag{3}$$

where the first condition Equation (2) is the normalizing condition and the second one Equation (3) is called the Studys' quadric [15]. The Study parameters may be obtained directly from a homogeneous transformation matrix using the following ratios [2–7]:

$$\begin{aligned} \mathbf{x}\_{0}: \mathbf{x}\_{1}: \mathbf{x}\_{2}: \mathbf{x}\_{3} &= 1 + a\_{11} + a\_{22} + a\_{33}: a\_{32} - a\_{23}: a\_{13} - a\_{31}: a\_{21} - a\_{12} \\ &= a\_{32} - a\_{23}: 1 + a\_{11} - a\_{22} - a\_{33}: a\_{12} + a\_{21}: a\_{31} + a\_{13} \\ &= a\_{13} - a\_{31}: a\_{12} + a\_{21}: 1 - a\_{11} + a\_{22} - a\_{33}: a\_{23} + a\_{32} \\ &= a\_{21} - a\_{12}: a\_{31} + a\_{13}: a\_{23} + a\_{32}: 1 - a\_{11} - a\_{22} + a\_{33} \end{aligned} \tag{4}$$

$$\begin{array}{l} y\_0 = -\frac{1}{2}(t\_z \mathbf{x}\_3 + t\_y \mathbf{x}\_2 + t\_x \mathbf{x}\_1) \\ y\_2 = -\frac{1}{2}(-t\_z \mathbf{x}\_1 - t\_y \mathbf{x}\_0 + t\_x \mathbf{x}\_3) \end{array} \tag{5}$$

$$\begin{array}{l} y\_2 = -\frac{1}{2}(-t\_z \mathbf{x}\_1 - t\_y \mathbf{x}\_0 + t\_x \mathbf{x}\_3) \\ y\_3 = -\frac{1}{2}(-t\_z \mathbf{x}\_0 + t\_y \mathbf{x}\_1 - t\_x \mathbf{x}\_2) \end{array} \tag{6}$$

where the terms *aij* are the entries of a 3 × 3 rotation matrix within a homogeneous 4 × 4 matrix and [*tx*, *ty*, *tz*] <sup>T</sup> is the translation vector from the homogeneous matrix. As seen in Equation (4) there are four ways to write the Study parameters and it is guaranteed that at least one way yields a nonsingular representation of the Euclidean displacement (i.e., the Study parameters are free of parametric singularities). As an example, as long as the Study parameters do not represent a rotation by a value of π, the first ratio from Equation (4) may be used [2].

In the context of mechanism analysis, the relative position between a mobile coordinate frame, which is located on the robot mobile platform (or end-effector) and a fixed coordinate frame (which is at the base of robot) is studied. Homogeneous transformation matrices may be used in this context (which are also called the Denavit-Hartenberg convention [16]) and the algorithm of applying this method is straightforward: the first step is to define each transformation matrix for each joint and link; the second step is to multiply the transformation matrices starting from the base coordinate frame towards the mobile coordinate frame (for robots that contain multiple kinematic chains this procedure is done for every chain). The Study parameters may be obtained thereafter from the transformation matrix using Equations (4) and (5). The trigonometric functions (sines and cosines) may be substituted with the tangent of half angle formulae resulting in a rational form of the Study parameters which may be simplified further (by factoring the denominator and dividing through the greatest common divisor as in [7]) to obtain polynomials (achieving thus an algebraic description of the constraint equations). A more direct method to obtain the constraint equations of a mechanism is to write directly the Study parameters for each element of the kinematic chain (joints and links) and then to use quaternion multiplication to obtain the constraints (this is analog to multiplying the homogeneous matrices). As mentioned before, the Study parameters are based on a dual quaternion which has a general form:

$$Q = \mathbf{x}\_0 + \mathbf{x}\_1 i + \mathbf{x}\_2 j + \mathbf{x}\_3 k + \frac{1}{2}\varepsilon (y\_0 + y\_1 i + y\_2 j + y\_3 k) \tag{6}$$

where *i*, *j*, *k* are the imaginary parts and the following multiplication rules apply: *i* <sup>2</sup> = *j* <sup>2</sup> = *k*<sup>2</sup> = <sup>−</sup>1, *ij* = *k*, *jk* = *i*, *ki* = *j*, *ji* = −*k*, *kj* = −*i*, *ik* = −*j*, ε - 0, ε<sup>2</sup> = 0. The conjugate of the dual quaternion (which represent the inverse transformation) is defined as [2,9]:

$$\overline{Q} = \mathbf{x}\_0 - \mathbf{x}\_1 i - \mathbf{x}\_2 j - \mathbf{x}\_3 k + \frac{1}{2} \varepsilon (y\_0 - y\_1 i - y\_2 j - y\_3 k); \qquad \mathbf{Q} \cdot \overline{\mathbf{Q}} = [1:0:0:0:0:0:0:0] \tag{7}$$

#### *The Study Parameters for MIS Procedures*

For the MIS procedures the two most fundamental instrument manipulations are: (i) the *linear insertion*/*retraction* of a medical tool (e.g., needle); (ii) the *Remote Center of Motion* (RCM) manipulation which is a combination of the insertion of a medical instrument with the orientation of the instrument where the entry-point remains unchanged. Consequently the insertion/retraction procedure has 1 DOF, whereas the RCM has a maximum of 4 DOF (3 orientations and one translation). The Study parameters

equations which describe these two procedures are based on dual quaternion multiplication and are presented further.

Let *P*(*x*0:*x*1:*x*2:*x*3:*y*0:*y*1:*y*2:*y*3) be a point on the Study quadric *x*0*y*<sup>0</sup> + *x*1*y*<sup>1</sup> + *x*2*y*<sup>2</sup> + *x*3*y*<sup>3</sup> = 0 which describe the position and orientation of the mobile platform (the end-effector) under the constraint of a robotic device:

(1) The *linear insertion*/*retraction* of a medical tool is achieved by the translation on a given axis and is defined by:

$$P\_{\rm ins} = P \cdot P\_{\rm l} \tag{8}$$

$$P = P\_{\rm inv} \cdot \overline{P}\_{l} \tag{9}$$

where the *Pins* represents the coordinates of the target point (encoding also the orientation of the medical tool) *P* represents the insertion point and *Pl* represents the insertion length; [1:0:0:0:0:−*l*/2:0:0] for the insertion along the *X*' axis (of the mobile frame), [1:0:0:0:0:0:−*l*/2:0] for the insertion along the *Y*' axis and [1:0:0:0:0:0:0:−*l*/2] for the insertion along the *Z*' axis. Equation (8) represents the linear insertion whereas Equation (9) represents the retraction (i.e., the inverse operation) and it is given by multiplying with the conjugate of the dual quaternion.

(2) The *RCM manipulation* between two configurations of the robot is given by:

$$P\_{\rm RCM} = P \cdot P\_{\rm o} \cdot P\_{\rm l} \tag{10}$$

$$P = P\_{\rm RCM} \cdot \vec{P}\_l \cdot \vec{P}\_o \tag{11}$$

where *PRCM* represents the Study parameters for a general point that is manipulated using the RCM concept, *Po* are the Study parameters which describe a change in orientation and *Pl* are the Study parameters describing the insertion of the medical instrument. Assuming that *u* is the tangent of the half angle of a rotation angle, the Study parameters [1:*u*:0:0:0:0:0:0], [1:0:*u*:0:0:0:0:0] and [1:*0*:0:*u*:0:0:0:0] represent rotations (by a value of *u*) about the *X*', *Y*' and *Z*' axes respectively. Equation (10) describes the manipulation using RCM whereas Equation (11) describes the inverse operation which is achieved by multiplying with the quaternion conjugate. To describe a general way in which the medical tool is displaced using RCM both Equations (10) and (11) are used:

$$P\_{R\gets M} = P \cdot P\_{o,1} \cdot P\_{l,1} \cdot P\_{l,1} \cdot P\_{o,1} \cdot P\_{o,2} \cdot P\_{l,2} \cdot P\_{l,2} \cdot P\_{o,2} \cdot \dots \cdot P\_{o,n} \cdot P\_{l,n} \tag{12}$$

Although the quaternions offer a nice analytic method to operate the constraint equations, the Euclidean representation of specific points (e.g., the insertion point) is also necessary to create an intuitive description of the robotic assisted medical act (i.e., the Cartesian coordinates are used for an unambiguous robot interface targeting medical experts, whereas the quaternions are used for the computation). The relations between the Study parameters and points in 3D Cartesian space may be obtained from Equation (5) together with the strong normalizing condition in Equation (1) and are [2–4]:

$$\begin{aligned} t\_x &= -2(x\_0y\_1 - x\_1y\_0 + x\_2y\_3 - x\_3y\_2) \\ t\_y &= -2(x\_0y\_2 - x\_1y\_3 - x\_2y\_0 + x\_3y\_1) \\ t\_z &= -2(x\_0y\_3 + x\_1y\_2 - x\_2y\_1 - x\_3y\_0) \end{aligned} \tag{13}$$

where again the *t* = [*tx, ty, tz*] <sup>T</sup> is the translation vector which describes the position of mobile frame origin relative to the fixed frame origin.

#### **3. Pro-Hep-LCT Parallel Robotic System**

The Pro-Hep-LCT parallel robotic system is designed for the minimally invasive procedure of the palliative treatment of the unresectable liver tumors using brachytherapy, chemotherapeutic agent delivery or ablation. The procedure requires the accurate insertion of brachytherapy needles on linear trajectories (and releasing them afterwards since the brachytherapy treatment is delivered in specialized rooms), under the real live imaging given by an intra-operatory ultrasound probe. Consequently, two independent modules are needed, one for the linear insertion of the needles, and the second one for the RCM manipulation of the intra-operatory ultrasound probe [17]. In previous work, the authors developed medical instruments with redundant DOFs for the needle insertion procedures (such as the ones in [18,19]). The use of such instruments increases the precision of insertion, and reduces the risk of the medical procedure. Such medical instruments will also be integrated on the Pro-Hep-LCT robotic system, namely: one brachytherapy instrument which will insert and release the brachytherapy needle with a redundant DOF; one US probe instrument which will insert/retract the probe (with a redundant DOF).

#### *3.1. Pro-Hep-LCT Parallel Robotic System Description*

The two guiding modules of the Pro-Hep-LCT robotic system are identical and operate "in mirror" relative to each other (as illustrated in Figure 1). This allows more access possibilities for the medical tools since both of the guiding modules may guide either the needle instrument or the ultrasound probe instrument.

**Figure 1.** The CAD model for the Pro-Hep-LCT robotic system.

The kinematic scheme of the intra-operatory US probe guiding module is presented in Figure 2. The robotic module contains two planar modules which guide the robot end effector through a pair of universal joints. Since the two planar mechanisms are displaced by a constant value, the end effector must have a free prismatic joint to allow its orientation. To facilitate the computation, the ultrasound probe guiding module is decoupled into 3 mechanisms and the relation among them is the following:


4. The last active DOF is actuated by the active joint *q*<sup>1</sup> and it displaces the entire mechanism in the *Z* direction of the fixed coordinate system by moving the mechanical frame on which the two planar mechanisms are assembled. This DOF will be used only to guide the medical tool at the patient entry point and the manipulation of the medical instrument thereafter is done by the robot end-effector using redundant DOFs.

**Figure 2.** The kinematic scheme of the intra-operatory US probe guiding module.

#### *3.2. Pro-Hep-LCT Kinematics*

The algebraic method used in the paper, for the Pro-Hep-LCT analysis has the potential to offer a generalized solution for the kinematics (both forward and inverse) and shows a complete description of the workspace (together with the singularities). Although at times the constraint equations become complicated the method exploits the fact that the constraints are algebraic varieties and Groebner bases may be computed to obtain "fundamental" polynomials which describe the constraint varieties. The constraint equations for the Pro-Hep-LCT parallel robotic system are obtained for all three component mechanisms and the forward and inverse kinematics are obtained using the following stepwise procedures:


Moreover, since the robotic system uses end-effectors for the insertion of the medical tools, the mathematical description for the medical procedure using the quaternion viewpoint is straightforward and also easy to implement in the robot control. Assuming that the robotic system is in a configuration where the medical tool tip is at the patient insertion point, the Study parameters of the mobile platform are known. The manipulation of the medial tool constrained by RCM is achieved using Equation (12) and the active joint parameters are obtained using the *inverse kinematics* algorithm.

*Observation 1*: Although the Study parameters show computational advantages, there is a fundamental problem that must be overcome. For the displacement of the ultrasound probe from an initial configuration to the next configuration by using Equation (12) the Cartesian coordinates of the entry point will not change. However, in-between the two configurations, there is no guarantee that the entry point will not move. A numerical solution for this drawback is to use smaller displacements for the tool manipulations and interpolation to determine the tool path.

#### 3.2.1. The Upper Planar Mechanism Kinematics

The *upper planar mechanism* (Figure 2) is described by the constraints of the point *Pu* (located at the intersection axes of *U*1). Using dual quaternion multiplication (described by Equation (6) and the rules of the imaginary parts multiplication) the constraint equations for the point *Pu* are:

$$P\_{\mathbf{u}}(\mathbf{x}\_{i},\mathbf{y}\_{i}) \to \begin{bmatrix} \mathbf{x}\_{\mathrm{u}1} \cdot \mathbf{P}\_{\mathbf{u}2} \cdot \mathbf{P}\_{\mathbf{u}3} \cdot \mathbf{P}\_{\mathbf{u}4} \end{bmatrix}$$

$$\begin{array}{c} P\_{\mathbf{u}}(\mathbf{x}\_{i},\mathbf{y}\_{i}) \to \begin{bmatrix} \mathbf{x}\_{\mathrm{0}} \cdot \mathbf{x}\_{1} \cdot \mathbf{x}\_{2} \cdot \mathbf{x}\_{3} \cdot \mathbf{y}\_{0} \cdot \mathbf{y}\_{1} \cdot \mathbf{y}\_{2} \cdot \mathbf{y}\_{3} \end{bmatrix}^{T} = \begin{bmatrix} 1 \cdot 0 \cdot 0 \cdot 0 \cdot 0 \cdot 0 \cdot \left. -\frac{1}{2} (l\_{q5} + l\_{q4}) \right| \cdot -\frac{1}{2} q\_{5} \circ 0 \end{bmatrix}^{T} \\\ l\_{q4} = \sqrt{l\_{d4}^{2} - (q\_{5} - q\_{4})^{2}} \end{array} \tag{14}$$

where the dual quaternions used in obtaining the *Pu* map are detailed in Table 1.


**Table 1.** The quaternions used in defining the constraints of the upper planar chain.

#### 3.2.2. The Lower Planar Mechanism Kinematics

For the *lower planar mechanism* (Figure 2) the Study parameters are written for the two kinematic chains (both RRR) which are coupled at the point *Pd* (located at the midpoint of the link *d*). Therefore the point *Pd* must simultaneously satisfy two equations:

 *Pd* <sup>=</sup> *Pd*<sup>1</sup> · *Pd*<sup>2</sup> · *Pd*<sup>3</sup> · *Pd*<sup>4</sup> · *Pd*<sup>5</sup> · *Pd*<sup>6</sup> *Pd* = *Pd*<sup>7</sup> · *Pd*<sup>8</sup> · *Pd*<sup>9</sup> · *Pd*<sup>10</sup> · *Pd*<sup>11</sup> · *Pd*<sup>12</sup> *Pd*(*xi*, *yi*) → [*x*<sup>0</sup> : *x*<sup>1</sup> : *x*<sup>2</sup> : *x*<sup>3</sup> : *y*<sup>0</sup> : *y*<sup>1</sup> : *y*<sup>2</sup> : *y*3] *<sup>T</sup>* = [*x*<sup>0</sup> :0:0: *x*<sup>3</sup> :0: *y*<sup>1</sup> : *y*<sup>2</sup> : 0] *T x*<sup>0</sup> = 4(*t*1*t*<sup>2</sup> + *t*1*t*<sup>3</sup> + *t*2*t*<sup>3</sup> − 1) *x*<sup>0</sup> = 4(*w*1*w*<sup>2</sup> + *w*1*w*<sup>3</sup> + *w*2*w*<sup>3</sup> − 1) *x*<sup>3</sup> = 4(*t*1*t*2*t*<sup>3</sup> + *t*<sup>1</sup> + *t*<sup>2</sup> − *t*3) *x*<sup>3</sup> = 4(*w*1*w*2*w*<sup>3</sup> + *w*<sup>1</sup> + *w*<sup>2</sup> − *w*3) *y*<sup>1</sup> = *t*1*t*2(*dt*<sup>3</sup> + 2*l*<sup>1</sup> − 2*l*2) + *t*1*t*3(2*l*<sup>1</sup> + 2*l*2) + *t*2*t*3(−2*l*<sup>1</sup> + 2*l*2) − *d*(*t*<sup>1</sup> + *t*<sup>2</sup> + *t*3) + 2(*l*<sup>1</sup> + *l*2) *y*<sup>1</sup> = *w*1*w*2*w*3(*d* + 2*lq*3) + 2*w*1*w*2(*l*<sup>3</sup> − *l*4) + 2*w*1*w*3(*l*<sup>3</sup> + *l*4) + 2*w*2*w*3(−*l*<sup>3</sup> + *l*4) − *d*(*w*<sup>1</sup> + *w*<sup>2</sup> + *w*3)+ 2*lq*4(*w*<sup>1</sup> + *w*<sup>2</sup> + *w*3) + 2(*l*<sup>3</sup> + *l*4) *y*<sup>2</sup> = 2*t*1*t*2*t*3(*l*<sup>2</sup> − *l*1) − *d*(*t*1*t*<sup>2</sup> + *t*2*t*<sup>3</sup> + *t*1*t*3) + 2*l*1(*t*<sup>1</sup> − *t*<sup>2</sup> − *t*3) + 2*l*2(*t*<sup>1</sup> + *t*<sup>2</sup> − *t*3) + *d y*<sup>1</sup> = *w*1*w*2*w*3(2*l*<sup>4</sup> − 2*l*3) − *w*1*w*2(*d* + 2*lq*4) − *w*1*w*3(*d* + 2*lq*4) − 2*w*2*w*3(*d* + 2*lq*4) + *l*3(*w*<sup>1</sup> − *w*<sup>2</sup> − *w*3)+ 2*l*4(*w*<sup>1</sup> + *w*<sup>2</sup> − *w*3) + *d* + 2*q*<sup>4</sup> (15)

The twelve quaternions used in obtaining the map *Pd* are detailed in Table 2. It is already visible from the zeroes of the Study parameters in Equation (15) that the point *Pd* is constrained to move on a plane. However, Equation (15) describe the constraints of a six bar mechanism and a change in orientation (of the link *d*) is still possible. The solution to this is to simply set *x*<sup>0</sup> = 1 and *x*<sup>3</sup> = 0, which will also satisfy the normalization condition defined in Equation (2).


**Table 2.** The quaternions used in defining the constraints of the lower planar chain.

<sup>1</sup> The *ti* and *wi* represent the tangent of half angles parameters.

The fact that the Study parameters Equation (15) describe varieties in the Study parameter space, allows the elimination of the *t*2, *t*3, *w*2, *w*<sup>3</sup> parameters (which are the free motion parameters in the planar mechanisms) by computing elimination ideals with Groebner bases (this technique is detailed in [20] and was used in the analysis of various robotic systems such as [4,7]); it is worth mentioning that the parameter elimination can be also achieved with the Linear Implicitization Algorithm (LIA) [21] which was used in [6]. The resulting constraints in this case will be the union of the equations found in the two Groebner bases which were computed to eliminate the free motion parameters and are:

$$\begin{array}{l} \mathbf{G}\_{d} = <16(l\_{1}^{2}+1)y\_{1}^{2}-16l\_{1}(l\_{1}^{2}-1)y\_{1}+16(l\_{1}^{2}+1)y\_{2}^{2}+8(dl\_{1}^{2}+4l\_{1}l\_{1}+d)y\_{2}+(d^{2}+4l\_{1}^{2}-4l\_{2}^{2})l\_{1}^{2}+8dl\_{1}l\_{1} \\ > & +d^{2}+4(l\_{1}^{2}-l\_{2}^{2});16(w\_{1}^{2}+1)y\_{1}^{2}-16l\_{3}(w\_{1}^{2}-1)y\_{1}+16(w\_{1}^{2}+1)y\_{2}^{2}+8((2l\_{4}d+d)w\_{1}^{2}+4l\_{3}w\_{1}+2l\_{4}d+d)y\_{2} \\ > & +(d^{2}+4(dl\_{q4}+l\_{3}^{2}-l\_{4}^{2}+l\_{q4}^{2}))w\_{1}^{2}+8l\_{3}(d+2l\_{4})w\_{1}+d^{2}+4(dl\_{q4}+l\_{3}^{2}-l\_{4}^{2}+l\_{q4}^{2})> \end{array} \tag{16}$$

#### 3.2.3. The Mobile Platform Mechanism Kinematics

For the *mobile platform* (Figure 2) the forward kinematics is written using the Study parameters with the observation that the inputs are the two points *Pu(Xu,Yu)* and *Pd(Xd,Yd)* (were the eight dual quaternions used in the *P* map are presented in Table 3):

$$\begin{aligned} \left\{ \begin{array}{l} P = P\_1 \cdot P\_2 \cdot P\_3 \cdot P\_4\\ P = P\_5 \cdot P\_6 \cdot P\_7 \cdot P\_8 \end{array} \right. \\ P\_u(\mathbf{x}\_i, y\_i) \to \left[ \mathbf{x}\_0 : \mathbf{x}\_1 : \mathbf{x}\_2 : \mathbf{x}\_3 : y\_0 : y\_1 : y\_2 : y\_3 \right]^T &= \left[ 2 : 2u\_1 : -2u\_2 : -2u\_1u\_2 : y\_0 : y\_1 : y\_2 : y\_3 \right]^T \\ \mathbf{y} = \left[ 2 : 2v\_1 : -2v\_2 : -2v\_1u\_2 : y\_0 : y\_1 : y\_2 : y\_3 \right]^T \\ y\_0 = -D\_n u\_1 u\_2 + X\_d u\_1 - Y\_d u\_2 + c u\_1 & y\_0 = v\_1 v\_2 (D\_n - D\_2 + p) + X\_d v\_1 - Y\_d u\_2 \\ y\_1 = Y\_d u\_1 u\_2 + D\_n u\_2 - X\_d & y\_1 = v\_2 (Y\_n u\_1 - D\_n - D\_2 - p) - X\_u \\ y\_2 = -(X\_d + c)u\_1 u\_2 + D\_n u\_1 - Y\_d & y\_2 = v\_1 (X\_u u\_2 - D\_n - D\_2 - p) - Y\_u \\ y\_3 = (X\_d + c)u\_2 + Y\_d u\_1 - D\_n & y\_3 = X\_d v\_2 + Y\_u u\_1 + D\_n - D\_2 + p \end{aligned} \tag{17}$$


**Table 3.** The quaternions used in defining the constraints of the robot mobile platform.

<sup>1</sup> The *ui* and *vi* represent the tangent of half angles parameters.

The Study parameters equations found in Equation (17) represent a variety in Study space. Like in the previous case the free motion parameters are eliminated from the equations (the rotation parameters *u*1, *u*2, *v*1, *v*<sup>2</sup> and the translation parameter *p*) by computing elimination ideals using Groebner bases. Consequently, the resulting input/output equations relate the terms *Xd, Xu, Yd, Yu* with the Study parameters (the mobile platform coordinates and orientation). The constraints of the mobile platform are given by the union of the Study parameter equations of the two elimination ideals (describing *P*—see Equation (17)) which are four linear equations (in the Study parameters), nine quadratic equations, one cubic equation, the Study quadric and the normalizing condition (due to the complexity of the equations they are shown only in general form—e.g., g1 is a linear equation in the unknowns *x*3, *y*0, *y*1, *y*2, g5 is an quadratic equation in the unknowns *y*0, *y*1, *y*2, *y*<sup>3</sup> and so on):

*G* =< *g*1(*x*3, *y*0, *y*1, *y*2), *g*2(*x*2, *y*0, *y*1, *y*3), *g*3(*x*1, *y*2, *y*0, *y*3), *g*4(*x*0, *y*1, *y*2, *y*3), *g*<sup>2</sup> <sup>5</sup>(*y*0, *y*1, *y*2, *y*3), *g*2 6..10(*x*0, *<sup>x</sup>*1, *<sup>x</sup>*2, *<sup>x</sup>*3, *<sup>y</sup>*0, *<sup>y</sup>*1, *<sup>y</sup>*2, *<sup>y</sup>*3), *<sup>g</sup>*<sup>2</sup> 11..12(*x*0, *<sup>x</sup>*1, *<sup>x</sup>*2, *<sup>x</sup>*3), *<sup>g</sup>*<sup>2</sup> <sup>13</sup>(*x*0, *<sup>x</sup>*1, *<sup>x</sup>*2, *<sup>x</sup>*3, *<sup>y</sup>*1), *<sup>x</sup>*<sup>2</sup> <sup>0</sup> <sup>+</sup> *<sup>x</sup>*<sup>2</sup> <sup>1</sup> <sup>+</sup> *<sup>x</sup>*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>x</sup>*<sup>2</sup> <sup>3</sup> − 1, *x*<sup>0</sup> *y*<sup>0</sup> + *x*<sup>1</sup> *y*<sup>1</sup> + *x*<sup>2</sup> *y*<sup>2</sup> + *x*<sup>3</sup> *y*3, *g*<sup>3</sup> <sup>16</sup>(*x*0, *x*1, *x*2, *x*3, *y*0, *y*1, *y*2, *y*3) > (18)

At this point it can be checked that the basis *G* has the Hilbert dimension 0 which corresponds to the forward kinematic solution being points (not necessarily unique since for one set of input parameters the mechanism may pose in different ways). Moreover, the basis may be solved for both the Study parameters (when the active joints are inputs) achieving the forward kinematics and for the active joints (when the Study parameters are inputs) achieving the inverse kinematics.

#### *3.3. Pro-Hep-LCT Singularities*

This section describes the Pro-Hep-LCT robotic system singularities which are determined from the constraint equations obtained in the previous section. Consequently, the singularities are analyzed for each component mechanism individually. This is possible due to the fact that the robot is a hybrid (and the end effector is guided by two planar mechanisms). The inputs of the lower planar mechanism (*t*<sup>1</sup> and *w*<sup>1</sup> being the tangents of the half angles of *q*<sup>2</sup> and *q*<sup>3</sup> respectively) do not influence the point *Pu* (of the upper mechanism) and the inputs of the upper chain (*q*<sup>4</sup> and *q*5) do not influence the point *Pd*. For the input singularities, while the robot is rigid, an infinitesimal change in the inputs will produce no change in the outputs, whereas for the output singularities, an infinitesimal change in the output is possible even though there is no change in the input. The singularities for the robotic system are determined using the rank of the input and output Jacobians for the three presented mechanisms. If the Jacobians of any of the mechanism are rank deficient a singularity occurs. The input and output Jacobians were determined for each component mechanism and are:

$$\begin{aligned} \left[j\_{i,1} \cdot \dot{q}\_u + j\_{o,1} \cdot \dot{P}\_u = 0 \quad \right] \dot{j}\_{i,1} &= \begin{bmatrix} \frac{\partial \left(G\_u[1], G\_u[2]\right)}{\partial \left(q\_4, q\_5\right)} \end{bmatrix}; \quad \begin{aligned} \left[j\_{o,1} = \begin{bmatrix} \frac{\partial \left(G\_u[1], G\_u[2]\right)}{\partial \left(y\_1, y\_2\right)} \end{bmatrix} \right] \dot{q}\_u &= \begin{bmatrix} \dot{q}\_4, \dot{q}\_5 \end{bmatrix}^T; \\ \dot{P}\_u &= \begin{bmatrix} \dot{y}\_1, \dot{y}\_2 \end{bmatrix}^T \end{aligned} \tag{19}$$

$$\begin{aligned} \left[\boldsymbol{J}\_{i,2} \cdot \dot{\boldsymbol{q}}\_d + \boldsymbol{J}\_{o,2} \cdot \dot{\boldsymbol{P}}\_d = 0 \quad \boldsymbol{J}\_{i,2} = \begin{bmatrix} \frac{\partial (\boldsymbol{G}\_u[1], \boldsymbol{G}\_u[2])}{\partial (\boldsymbol{q}\_d \boldsymbol{q}\_5)} \end{bmatrix} \right] \quad \boldsymbol{J}\_{o,2} = \begin{bmatrix} \frac{\partial (\boldsymbol{G}\_d[1], \boldsymbol{G}\_d[2])}{\partial (\boldsymbol{y}\_1 \boldsymbol{y}\_2)} \end{bmatrix} \dot{\boldsymbol{q}}\_d = \begin{bmatrix} \dot{\boldsymbol{I}}\_1, \dot{\boldsymbol{w}}\_1 \end{bmatrix}^T \end{aligned} \tag{20}$$

$$\begin{aligned} \left[J\_{i,3} \cdot \dot{q} + J\_{o,3} \cdot \dot{P} = 0 \quad J\_{i,3} = \begin{bmatrix} \frac{\partial (G[1] \dots G[16])}{\partial (q\_4 \rho\_5)} \end{bmatrix} \right] \quad \begin{aligned} \left[J\_{o,3} = \begin{bmatrix} \frac{\partial (G^\*[1] \dots G^\*[8])}{\partial (\dot{x}\_i y\_i)} \end{bmatrix} \dot{q} = \begin{bmatrix} \dot{X}\_{\mu\_1} \dot{Y}\_{\mu\_1} \dot{X}\_{d\mu} \dot{Y}\_d \end{bmatrix}^T \\\ \dot{P}\_d = \begin{bmatrix} \dot{x}\_0 \dots \dot{x}\_3 \dot{y}\_0 \dots \dot{y}\_3 \end{bmatrix}^T \end{aligned} \tag{21}$$

The *G*\* basis represents a pure lexicographic Groebner base computed from the *G* base and has the advantage that it has 8 equations and the Jacobian *Jo*,3 remains a square matrix. This is not true for the *Ji*,3, therefore, all the determinants of the square minors must be computed for the singularity analysis.

Analyzing the rank deficiency conditions of the Jacobians, the output singularities are given by the following 12 polynomial equations (Equations (24)–(27) and (32)–(36) are only presented in the general form since they are lengthy):

$$\text{sg}\_1: w\_1^2 + 1 = 0 \tag{22}$$

$$\text{sg}\_2: t\_1^2 + 1 = 0 \tag{23}$$

$$\text{sg3}(l\_1, w\_1, l\_1, l\_3) = 0 \tag{24}$$

$$\text{sg}\_4(l\_{1\prime}w\_{1\prime}l\_{1\prime}l\_{3\prime}l\_{q4\prime}d) = 0\tag{25}$$

$$\log\_5(t\_{1\prime}w\_{1\prime}l\_{1\prime}l\_{2\prime}l\_{3\prime}l\_{q4\prime}d) = 0\tag{26}$$

$$\text{sg}\_6(l\_1, w\_1, l\_1, l\_2, l\_3, l\_{q4}, d) = 0 \tag{27}$$

$$\{\mathbf{s}\mathbf{g}\_7 \colon \mathbf{Y}d - \mathbf{Y}u = 0\tag{28}$$

$$s\!\!g\_8: Xd - Xu + c = 0\tag{29}$$

$$s\emptyset \!\!g : D\_{\mathbb{Z}} = 0$$

$$s g\_{10}: D\_z^2 + \chi\_d^2 - 2\Upsilon\_d \Upsilon\_u + \chi\_u^2 = 0 \tag{31}$$

$$s\!g\_{11}(X\_{\mu}, Y\_{\mu}, X\_{d\nu} \; Y\_{d\nu} D\_{\overline{z}}, \mathfrak{c}) = 0\tag{32}$$

$$\log\_{12}(X\_{\mu\prime}Y\_{\mu\prime}X\_{d\prime}Y\_{d\prime}D\_{z\prime}\mathfrak{c}) = 0\tag{33}$$

and the input singularities are given by three polynomials and an irrational function (due to the substitution found in the map defined in Equation (14)):

$$s\_{\mathbb{S}\backslash 3}(y\_1, y\_2, l\_{\mathbb{S}\backslash}l\_4, l\_{\mathbb{q}\backslash}d) = 0\tag{34}$$

$$s\_{\ $}g\_{14}(y\_1, y\_2, l\_{\$ }, l\_4, l\_{\mathfrak{q}4}, d) = 0\tag{35}$$

$$s\_{\mathbb{S}^1 5}(y\_1, y\_2, l\_{3\prime}, l\_{4\prime}, l\_{q4\prime}, d) = 0\tag{36}$$

$$\text{sg}\_{16} \colon \frac{q\_4 - q\_5}{\sqrt{(q\_4 - q\_5 + l\_{d4})(-q\_4 + q\_5 + l\_{d4})}} = 0 \tag{37}$$

Equations (22)–(37) represent the general singularities of the Pro-Hep-LCT robotic system. The simple equations are straight forward to analyze: Equations (22), (23) and (31) have no real solutions and are not of interest in kinematics; Equation (30) cannot be true since it implies that the two planar mechanisms operate in the same plane (which is not the case); Equation (37) represents a

singularity when active joints *q*<sup>4</sup> and *q*<sup>5</sup> are overlapped or are at a distance *lq*<sup>4</sup> apart (both conditions being physically impossible due to mechanism design considerations). Although the conditions from Equations (24), (28) and (29) make the Jacobians rank to drop, they are not singularities of the constraint varieties. Further investigations showed that the algorithm used in computing the Groebner bases *Gd* and *G* divided through (*t*<sup>1</sup> ± *w*1), (*Xd* − *Xu* + *c*) and (*Yd* − *Yu*) respectively. Consequently the bases are not correct when these conditions are met since the division by 0 is not allowed. In a sense, Equations (24), (28) and (29) represent singularities of the computed Groebner bases and are not singularities of the constraint varieties. To check that this is true, the conditions were substituted into the kinematic model before the computation of the Groebner bases (*Gd* and *G*) and no singularity was found.

For the lengthy polynomial equations which describe singularities of the Pro-Hep-LCT the numerical values for the geometric parameters resulted from the robot design (the links lengths as shown in Table 4) were used to evaluate the equations to allow a geometric interpretation of the singularities.

**Table 4.** Numerical values for the lengths of the Pro-Hep-LCT robotic system.


After the numerical evaluation of the remaining equations, Equations (25) and (27) showed only complex solutions and are not of interest. Equation (26) showed a configuration where the link *d* is aligned with either of the links *l*<sup>2</sup> or *l*<sup>4</sup> (see Figure 3a). This singularity represents the separation of two working modes of the lower planar mechanism and it is part of the boundary of the operational workspace (as illustrated in following sections). Equations (32) and (33) showed a configuration which was already described in Equation (30), namely the *Dz* parameter being zero. However, this condition is not allowed since the mechanism in not constructed in such way.

**Figure 3.** Singularities of the Pro-Hep-LCT robotic system; (**a**) configuration where the links l2 and d are aligned; (**b**) self-motion when the links l2 and l4 are parallel; (**c**) configuration were the links l1 and l2 are aligned.

For the input singularities the equations were also evaluated with the numerical values detailed in Table 4. Through a close analysis it was found that Equation (34) describes a self-motion of the link d which is present when the links *l*<sup>2</sup> and *l*<sup>4</sup> are parallel (Figure 3b). This configuration must be avoided using the robot control since it may cause harm to the patient during the medical procedure. Equation (35) describes a singularity when *t*<sup>1</sup> = 1 and *w*<sup>1</sup> = –1 which is physically impossible since in the robot design it implies that various mechanical links (*l*3, *l*4) will overlap. Equation (36) describes configurations where the link *l*<sup>1</sup> is aligned with *l*<sup>2</sup> (Figure 3c) or the link *l*<sup>3</sup> is aligned with *l*4. This actually represents the boundary of the workspace of the lower planar mechanism.

#### *3.4. Pro-Hep-LCT Workspace*

Analyzing the constraint equations for the Pro-Hep-LCT robotic system (Equations (14), (16) and (18)) it can be shown that the lower planar mechanism has four working modes, the upper planar mechanism has two working modes and even the end effector has four working modes (at least theoretically). For this reason, the geometry of the workspace is complex and the focus hereafter switches towards the operational workspace (the set of configurations intended for achieving the medical task).

The operational workspace was analyzed in two ways: at first the operational workspace describing the tip of the medical tool was determined (which for the medical task at hand represents the skin entry-point); and second, the redundant DOF for the medical tool was introduced to geometrically illustrate the Pro-Hep-LCT robotic system capabilities to achieve the medical task. Figure 4a shows the boundary of the operational workspace in which the end-effector of the Pro-Hep-LCT robotic system has no orientation (i.e., the medical tool points vertically downwards) and Figure 4b shows a cross-section of the space. Within this cross-section of the workspace, there exists a subspace (on the right hand side of the blue line *X* = 214, see Figure 3b) in which all the possible orientations of the medical tool are possible (i.e., the dexterous workspace). In this region of the workspace the robotic assisted medical procedure should be the least limited (due to the mechanism constraints). However, there is a point that describes the self-motion of the mechanism (of coordinates *X* = 286.66 mm, *Y* = 200 mm) which must be avoided (using the control algorithm) since in this configuration the control of the robot is lost and injury may occur.

**Figure 4.** Pro-Hep-LCT operational workspace; (**a**) envelope; (**b**) cross-section of the dexterous work-space.

To illustrate the operational workspace while the medical tool is inserted in the patient, two entry points are chosen. First, an entry point was chosen within the dexterous workspace (*X* = 260 mm, *Y* = 180 mm, *Z* = 200 mm) and Figure 5a illustrates two spherical surfaces defining two different depths of insertion (108 mm and 70 mm, respectively). The second example was chosen outside the dexterous workspace (see Figure 5b), having the entry point at (*X* = 200 mm, *Y* = 180 mm, *Z* = 200 mm), with two insertion depths (108 mm and 70 mm, respectively).

**Figure 5.** Pro-Hep-LCT operational workspace for the RCM manipulation; (**a**) entry point [260 mm, 180 mm, 200 mm]; (**b**) entry point [200 mm, 180 mm, 200 mm].

#### **4. Numerical Results and Discussion**

A numerical simulation is illustrated in this section to show how the results presented in Section 2 apply to the Pro-Hep-LCT robotic system functionality. To illustrate an example a Cartesian point is chosen as an entry-point *E* on the patient abdomen with the numerical coordinates [*XE* = 300 mm, *YE* = 180 mm, *ZE* = 110]. Furthermore, for this example, the ultrasound probe at the entry-point is chosen in a configuration where it points directly downwards in a vertical pose. The configuration has the following numerical values of the Study parameters [*x*<sup>0</sup> = 1, *x*<sup>1</sup> = 0, *x*<sup>2</sup> = 0, *x*<sup>3</sup> = 0, *y*<sup>0</sup> = 0, *y*<sup>1</sup> = −150, *y*<sup>2</sup> = −90, *y*<sup>3</sup> = 55] (relative to the fixed coordinate frame *XYZ;* to achieve this *q*<sup>1</sup> = 220 mm is a necessary condition). The inverse kinematics is computed from Equations (14), (16) and (18) with the numerical values for the Study parameters substituted in, yielding the numerical values for the active joint parameters [*t*<sup>1</sup> = −0.1343, *w*<sup>1</sup> = 0.0214, *q*<sup>4</sup> = 20.802 mm, *q*<sup>5</sup> = 180 mm] (solution describing the intended working mode for the Pro-Hep-LCT robotic system). Table 5 shows the initial configuration of the ultrasound probe (O at the entry-point) and two points A and B which are the end points of two simple tool paths (at 100 mm insertion depths). The Study parameters required for the RCM manipulation are also shown for the two tool paths with a discretization by angular values of 5◦ starting from 0◦ (the initial configuration) and ending at 30◦ (for point *A* with ϕ being a rotation about *X'* axis and for point *B* where θ represents a rotation about *Y*' axis).


**Table 5.** Numerical values for two paths of the ultrasound probe constrained by RCM.

The *t*<sup>1</sup> and *w*<sup>1</sup> are the tangent of the half angle of *q*<sup>2</sup> and *q*<sup>3</sup> respectively; n represents the value for the redundant DOF from the robot end-effector.

Lastly, Table 5 shows the numerical values of the active joint parameters for each ultrasound probe pose (computed by means of inverse kinematics using Equations (14), (16) and (18)).

Figure 6 illustrates the discretization of the two paths of the ultrasound probe tip starting from a point O' (where the tool is in a vertical pose and inserted at a 100 mm depth) and ending at the points A and B respectively (with 30◦ angle values about *X*' and *Y*' axes respectively). These paths are simulated with the ultrasound probe constrained by the RCM. Figure 6a illustrates the Cartesian positions of each computed point from the two paths whereas Figure 6b illustrates the active joint values changing over a period of 6 s (using linear interpolation) required for the specified paths.

**Figure 6.** Ultrasound probe tip path simulation: (**a**) path discretization; (**b**) active joint space.

The computed numerical values (for the active joints) that describe the two paths were used in a kinematic simulation using the Siemens NX software. The simulation purpose was to analyze the motion of the insertion point in Cartesian coordinates during the ultrasound probe manipulation on the two paths *A* and *B* over a period of 6 s (since in between two consecutive points of a path there is no guarantee that the RCM point is unchanged). For the path O'A, Figure 7a shows the simulation from the Siemens NX environment where the variation of the entry-point E in each *X*, *Y* and *Z* direction relative to the predefined point [*XE* = 300 mm, *YE* = 180 mm, *ZE* = 100] is highlighted. On the *X* direction, the maximum variation during the simulation was about 0.5 mm, reaching to a maximum error of 1.45 mm. (error that may be introduced by numerical errors of the interpolation). On the *Y* direction, the maximum variation was 0.05 mm adding to a maximum error of 0.15 mm which is less than the total error from the *X* displacement.

On the Z direction the maximum variation was about 0.12 mm. leading to a total error of 0.18 mm. For the path O'B, the simulation is shown in Figure 7b where the error in the displacement on *X* direction is about 3 mm, the error in the Y direction is about 2 mm and the error in the *Z* direction is about 0.2 mm. There are two types of fluctuations present in the displacement graphs illustrated in Figure 7, the first type being most likely due to the kinematic model and it was expected (see Observation 1 in Section 3.1) since the RCM is guaranteed to be kept only in the control points of the path and not in between. The second fluctuation refers to the shift in the curve (e.g., the 5 s mark for the *Y* displacement in Figure 7a) which is introduced by interpolation errors, most likely amplified by changing the simulation environment from Maple (where the inverse kinematics was computed) to Siemens NN (where the RCM point was studied). However, the magnitude of the errors presented so far is not significant for the medical task since a bigger error may be introduced by the tissue motion during the respiration process.

**Figure 7.** Simulation in the Siemens NX environment (variation of the entry-point Cartesian coordinates); (**a**) displacement of the entry point while the ultrasound probe was manipulated on path O'A; (**b**) displacement of the entry point while the ultrasound probe was manipulated on path O'B.

By utilizing the Study parameters method the authors managed to achieve a complete geometric analysis (obtaining the constraints, singularities and describing the workspace) of the Pro-Hep-LCT parallel robotic system designed for the minimally invasive task of targeted liver treatment (through brachytherapy, ablation of chemotherapy) under ultra sound imaging. Although the robotic system may appear complex, analyzing the "segments" that compose it (referred as the lower, upper and mobile platform mechanisms in Section 3) facilitated the computation of the constraints and singularities without any disadvantages (in the beginning of Section 3.2 it is explained why segmenting the robotic system is allowed). Moreover, the use of quaternion multiplication for the medical tools manipulation (Equations (8)–(12)) offered a distinct computation advantage since it complements well with the constraints computed with the Study parameters model. The dual quaternion approach for RCM manipulation is not novel (see e.g., [9]) and it has been shown that the main disadvantage (which is discussed in Observation 1 Section 3.1) may be minimized using interpolation (this method was proven on a serial manipulator in [9]). Since the solutions for the forward kinematics of the Pro-Hep-LCT parallel robotic system had the dual quaternion form (i.e., Study parameters) the manipulation of the medical instruments constrained by RCM was straightforward (see Section 2 Equations (8)–(12)). By means of linear interpolation it was also shown that the displacement of the insertion point is insignificant for the medical procedure. Of course more refined interpolation methods (or more control points) may be used if needed to increase precision but the tradeoff must be considered (computation time vs. precision since a real time control must be achieved). Besides RCM manipulation using dual quaternions, other medical applications may be considered as well, e.g., motor rehabilitation, or even determining the kinematics of the human finger in an analytical manner (in contrast to the experimental models—see e.g., [22]).

#### **5. Conclusions**

Based on a mathematical method that describes the global kinematics (using the Study parameters of the total Euclidean displacement), the Pro-Hep-LCT robotic system was completely analyzed to determine its forward and inverse kinematics, singularities and workspace. To determine the robot capabilities to perform the medical task, the operational workspace of the robotic system was described in two different ways: first, based on the medical relevant position of the robot (relative to the patient) the singularity free operational workspace was shown which also contains the dexterous workspace (which is especially relevant for the medical task); second, the manipulation of the medical instrument under the constraint of the RCM was studied creating also a description of the active joint space in a discrete manner. The analysis of the singularity free operational workspace shows that the Pro-Hep-LCT parallel robotic system is feasible for the RCM manipulation within the angular values of ±30◦ (relative to the vertical direction) in all directions. The numerical values of the active joint parameters which were computed via inverse kinematics (for consecutive points on predefined tool paths) were linearly interpolated and used as inputs for a kinematic simulation in the Siemens NX environment which showed a maximum relative deviation (from the insertion point) of 1.45 mm. The preliminary numerical results show that the robot control based on the Study parameters is feasible due to the insignificant deviation of the insertion point (with respect to the imposed point).

Further work is intended for developing the robotic system, namely the development of the mechanical architecture of the robotic system, the development of a controls system which uses the constraint equations together with the instrument manipulation method based on the dual quaternions. Moreover, further work is intended to test the precision of the prototype and optimize the control to account for needed corrections during the RCM manipulation (e.g., when parasite motions are present).

**Author Contributions:** Investigation, I.B., M.H., C.V., N.P., A.N. and D.P.; Methodology, I.B., M.H., C.V., N.P., A.N. and D.P.; Project administration, D.P.; Validation, M.H. and N.P.; Writing—original draft, I.B.; Writing—review & editing, I.B., C.V. and D.P.

**Funding:** This work was supported by a grant of the Romanian Ministry of Research and Innovation, PCCCDI – UEFISCDI, project number PN-III-P1-1.2-PCCDI-2017-0221/59PCCDI/2018 (IMPROVE), within PNCDI III and by project ExNanoMat, contract no. 21 PFE-2018 within PNCDI III.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
