Design and Visual Servoing Control of a Hybrid Dual-Segment Flexible Neurosurgical Robot for Intraventricular Biopsy

The MicroNeuro is designed for intraventricular neurosurgery. Based on knowledge of brain anatomy and clinical demand from surgeons, the main design goals are first summarized as follows:

1. 1.

   Dimension: The mean diameters of the foramen of Monro (FM) were 5.7 mm on the axial image, 7.8 mm on the coronal image, and 5.6 mm on the sagittal image [22]. Thus, the outer diameter of the flexible endoscope should be less than 5.4 mm to avoid collision with the FM.

2. 2.

   Endoscope features: The MicroNeuro should provide high quality images and a working channel for biopsy instruments. Since clinical surgery needs to be performed underwater, the MicroNeuro also needs to provide irrigation and suction functions.

3. 3.

   Dexterity: Deflective length of the MicroNeuro should be short and able to bend with a large curvature.

This work was developed based on the surgical robot system for neurosurgery, designated MicroNeuro [23]. As shown in Fig. 2(a), this system mainly consists of the MicroNeuro and its actuation units, which are mounted on the end of a 7 DoFs robot arm (ER7 Pro, ROKAE). The quick release mechanism of the MicroNeuro facilitates the individual disinfection of endoscopes. Besides, a control console is also built for master-slave teleoperation with four monitors, a foot pedal, a joystick (TCA, THRUSTMASTER) and a master device (TouchX, 3D SYSTEM).

The MicroNeuro consists of two bendable flexible robots which are connected to a rigid tube. As shown in Fig. 2(d), (e) and (f), it provides several functions, such as multi-view images, water irrigation and suction, working channel (diameter 1.2mm) and illumination. The distal end of the inner endoscope and the rigid catheter are each equipped with a camera (OV6946). Unlike conventional dual-segment flexible robots with fixed length, each robot of the MicroNeuro can be axially translated relative to each other, so two combined bending modes can be realized: (i) mode 1 [see Fig. 2(b)], the inner endoscope has no axial movement, and only the outer flexible sheath bends; (ii) mode 2 [see Fig. 2(c)], the inner endoscope could be inserted independently (maximum distance is 40mm).

The backbones of each flexible robot are manufactured by femtosecond laser cutting of superelastic nitinol tubes. Fig. 3 shows the parameter definitions and values. The two robots have multiple pairs of notched joints distributed along the axial direction, and each joint has a bidirectional symmetrical rectangular notch. Three nitinol cables, driven by brushless coreless motors (ASSUN), are welded to the distal end of each flexible robot and routed along a crimped grooves.

The distribution of notches in the backbone makes it axial stiffness larger than that in lateral direction, so the backbone would bend when the eccentrically fixed cables are stretched. Referring to the piecewise constant curvature (PCC) model [24], each segment of MicroNeuro bends with a constant curvature along its length, similar to a circular arc, when actuated. As shown in Fig. 3, MicroNeuro can be geometrically parameterized by $\bm{{\Phi}}=\left(\begin{matrix}z_{b}&\theta_{s}&\varphi_{s}&z_{e}&\theta_{e}&\varphi_{e}\end{matrix}\right)^{\mathrm{T}}$ in the configuration space, where ${z_{b}}$ is the overall insertion distance provided by the robot arm, $\theta_{s}$ and $\theta_{e}$ are the bending angles, $\varphi_{s}$ and $\varphi_{e}$ are the rotation angles between the bending plane and the $oxz$ plane, and $z_{e}$ is the variable length of the inner endoscope, provided by the servo motors. $\theta_{s}$, $\varphi_{s}$, $\theta_{e}$, $\varphi_{e}$ can be calculated from the actuator space variables $\bf q=\left(\begin{matrix}z_{b}&l_{s,1}&l_{s,2}&l_{s,3}&z_{e}&l_{e,1}&l_{e,2}&l_{e,3}\end{matrix}\right)^{\mathrm{T}}$:

$$\begin{split}{\theta_{i}}&={\frac{2\sqrt{l^{2}_{i,1}+l^{2}_{i,2}+l^{2}_{i,3}-l_{i1}l_{i2}-l_{i2}l_{i3}-l_{i1}l_{i3}}}{3\rho_{i}}}\\ {\varphi_{i}}&=\mathrm{tan}2({{l_{i1}+l_{i3}-2l_{i2}},{\sqrt{3}(l_{i3}-l_{i1})}})\end{split}$$

(1)

where $i\in\left\{e,s\right\}$, and the subscripts $e$ and $s$ used to represent the outer sheath and inner endoscope, respectively, $\rho_{i}$ is the distance between the center of the cable and the center of the robot, $l_{i,m}$, $m\in\left\{1,2,3\right\}$ are the length of the driving guide wires in each flexible robot. The transformation matrix ${{}^{b}_{{et}}}{{\mathrm{\bf T}}}\in{\mathbb{R}}^{4\times 4}$ from the base frame $O_{b}X_{b}Y_{b}Z_{b}$ to the robot tip frame $O_{et}X_{et}Y_{et}Z_{et}$ is:

$$\begin{split}{}_{{\bf{et}}}^{\bf{b}}{\bf{T}}=&{{\bf{\tau}}_{\bf{z}}}({z_{b}}){{\bf{R}}_{\bf{z}}}({\varphi_{s}}){{\bf{\tau}}_{\bf{x}}}(\frac{{{z_{s}}}}{{{\theta_{s}}}}){{\bf{R}}_{\bf{y}}}({\theta_{s}}){{\bf{\tau}}_{\bf{x}}}(-\frac{{{z_{s}}}}{{{\theta_{s}}}})\\ &{{\bf{R}}_{\bf{z}}}({\varphi_{e}}){{\bf{\tau}}_{\bf{x}}}(\frac{{{z_{e}}}}{{{\theta_{e}}}}){{\bf{R}}_{\bf{y}}}({\theta_{e}}){{\bf{\tau}}_{\bf{x}}}(-\frac{{{z_{e}}}}{{{\theta_{e}}}})\end{split}$$

(2)

where ${{\bf{\tau}}_{\bf{j}}}$, ${\bf R_{j}}\in{\mathbb{R}}^{4\times 4}$ respectively denote translation and rotation about axis j, $z_{s}$ is the length of the outer robot. Considering the offset $d$ of the camera frame $O_{c}X_{c}Y_{c}Z_{c}$ from the robot tip, the camera w.r.t. the base is

$$\begin{split}{}_{{\bf{c}}}^{\bf{b}}{\bf{T}}={}_{{\bf{et}}}^{\bf{b}}{\bf{T}}{{\bf{R}}_{\bf{z}}}(-{\varphi_{s}}-{\varphi_{e}}){{\bf{\tau}}_{\bf{y}}}(-d)\end{split}$$

(3)

The Jacobian matrix ${{\bf{J}}_{\bf{r}}}\in{\mathbb{R}}^{3\times 6}$ is used to analytically establish the approximate relationship between camera velocity and joint velocity. Considering the translational motion, at discrete instance $k$, the iterative form is $\Delta{\bf{P_{r}}}(k)={\bf{J_{r}}}(k)\Delta{\bf{\Phi}}(k)$, where $\Delta{\bf{P}}(k)={\bf{P}}(k+1)-{\bf{P}}(k)$ is the small displacement of the camera, ${\bf{\Phi}}(k)={\bf{\Phi}}(k+1)-{\bf{\Phi}}(k)$ and ${\bf{J_{r}}}(k)$ can be derived through forward kinematics ${}_{{\bf{et}}}^{\bf{b}}{\bf{T}}$. To reach a given target position ${{\bf{P}}_{\bf{G}}}\in{\mathbb{R}}^{3}$ of the end of the robot in $O_{et}X_{et}Y_{et}Z_{et}$, we need to inversely solve the appropriate joint configuration. The damped least squares method [25] provides an alternative Jacobian matrix to avoid joint velocity near singularities, i.e.

$$\begin{split}\Delta{\bf{\Phi}}(k)={\bf{J^{T}_{r}}}(k)({\bf{J_{r}}}(k){\bf{J^{T}_{r}}}(k)+\sigma{I})^{-1}({\bf{P_{G}}}(k)-{\bf{P}}(k))\end{split}$$

(4)

However, material nonlinearity, segment interaction, external loads, etc. may have a significant negative impact on the accuracy of the PCC model. In this work, we consider a moving camera while the targets are fixed at any instance $k$. As shown in Fig. 3(e), for a given point ${\bf{A}}\in\mathbb{R}^{3}$ in $O_{c}X_{c}Y_{c}Z_{c}$, its coordinates in the image frame $O_{I}xy$ and pixel frame $O_{p}uv$ are ${\bf{A}}(k)=(x,y)^{\mathrm{T}}$ and ${\bf{\varsigma}}(k)=(u,v)^{\mathrm{T}}$, respectively. According to the pinhole camera model, the perspective equation can be obtained from the relationship on similar triangles, i.e.

$$\begin{split}u=\frac{\lambda_{x}x}{\lambda}+c_{c},\enskip v=\frac{\lambda_{y}y}{\lambda}+c_{y}\end{split}$$

(5)

The motion of the features $\Delta\bm{\varsigma}(k)$ on the pixel plane can be predicted using the interaction matrix:

$$\begin{split}\Delta\bm{\varsigma}(k)={\bf{L_{m}}}(k)\Delta{\bf{P}}(k)\end{split}$$

(6)

where ${\bf{L_{m}}}\in\mathbb{R}^{2\times 3}$ is a block matrix of ${\bf{L_{o}}}=[{\bf{L_{m}}}^{2\times 3}|{\bf{L_{\omega}}}^{2\times 3}]$ related to linear velocity, and

$${\bf{L_{o}}}=\begin{bmatrix}-\frac{\lambda}{z_{c}}&0&\frac{x}{z_{c}}&\frac{xy}{\lambda}&-\frac{\lambda^{2}+x^{2}}{\lambda}&y\\ 0&-\frac{\lambda}{z_{c}}&\frac{y}{z_{c}}&-\frac{\lambda^{2}+y^{2}}{\lambda}&-\frac{xy}{\lambda}&-x\end{bmatrix}$$

(7)

where $\lambda_{x}$, $\lambda_{y}$ are the focal length in pixels, $c_{x}$, $c_{y}$ are optical center in pixels and $\lambda$ is focal length in millimeter. Define ${\bf{J_{a}}}(k)\in\mathbb{R}^{6\times 8}$ as the Jacobian matrix between the actuator space and the configuration space from Eq. (1), that is, $\Delta{\bf{\Phi}}(k)={\bf{J_{a}}}(k)\Delta{\bf{q}}(k)$. Combined Eq. (4) and (6), the overall Jacobian matrix between pixel velocity and actuator velocity can be derived as follow:

$$\Delta\bm{\varsigma}(k)={\bf{L_{m}}}(k){\bf{J_{r}}}(k){\bf{J_{a}}}(k)\Delta{\bf{q}}(k)$$

(8)

In the classic IBVS [26], there are several choices for the depth $z_{c}$ in the matrix ${\bf{L_{m}}}(k)$. In this study, the depth $z_{c}^{*}$ at the desired position was used, and ${\bf{\hat{L_{m}}}}(k)$ denotes the estimation matrix.

As a continuum robot, MicroNeuro has infinite DoFs. When subject to model mismatch problems caused by disturbance or manufacturing error, the model-dependent robot Jacobian matrix ${\bf{J}}(k)={\bf{J_{r}}}(k){\bf{J_{a}}}(k)$ may cause control deviations and need to be estimated online. First, the Jacobian estimate at $k=0$ needs to be obtained offline, then the Jacobian can be updated iteratively online during the robot movement.

1. 1.

   Initialization: A small actuator movement $\Delta{\bf{q_{+}}}(0)$ is imposed on the MicroNeuro while it is located outside the brain, and an external electromagnetic sensor (NDI Aurora) is mounted on the tip of MicroNeuro to measure the displacement. The $i$-th independent actuator variables $\Delta{{q_{+,i}}}(0)$ causes a position deviation of the camera $\Delta{\bf{P_{c,i}}}(0)$. Hence, ${\bf{\hat{J_{+}}}}(0)$ is constructed as:

   $${\bf{\hat{J_{+}}}}(0)=\left[\begin{matrix}\frac{\Delta{\bf{P_{c},0}}(0)}{\Delta{{q_{+,0}}}(0)}&\cdots&\frac{\Delta{\bf{P_{c},8}}(0)}{\Delta{{q_{+,8}}}(0)}\end{matrix}\right]$$

   (9)

   To reduce manufacturing error, ${\bf{\hat{J_{-}}}}(0)$ is similarly constructed while a opposite displacement $\Delta{\bf{q_{-}}}(0)=-\Delta{\bf{q_{+}}}(0)$ is imposed. ${\bf{\hat{J}}}(0)$ is set as:

   $${\bf{\hat{J}}}(0)=0.5({\bf{\hat{J_{+}}}}(0)+{\bf{\hat{J_{-}}}}(0))$$

   (10)

2. 2.

   Online Estimation: The alterations in the MicroNeuro position and Jacobian matrix between adjacent instance are small, thus, the current analytical Jacobian matrix ${\bf{{J}}}(k)$ could be appropriately adjusted using ${\bf{\hat{J}}}(k-1)$:

   $${\bf{\hat{J}}}(k)=(1-\omega(k)){\bf{J}}(k)+\omega(k){\bf{\hat{J}}}(k-1)$$

   (11)

   where $\omega(k)=\frac{1}{1+\epsilon(k)}$ is the weighting factor, and $\epsilon(k)=||\bm{\varsigma}(k)-{{\bm{\varsigma}}_{\bf{G}}}(k)||_{2}$ denotes as the distance between measured feature $\bm{\varsigma}(k)$ and the target feature ${{\bm{\varsigma}}_{\bf{G}}}(k)$. The normalized $\bm{\varsigma}(k)$ and ${{\bm{\varsigma}}_{\bf{G}}}(k)$ could be applied in $\omega(k)$.

The goal of the IBVS task is to minimize the error $\omega(k)$. Inspired from [27, 28], to reduce the negative impact of model inaccuracy and external disturbance, an internal model control (IMC) scheme [29] is applied in the visual MPC controller, as shown in Fig. 4. ${\bf{e}}(k)$ is defined as the predictive error, that is, ${\bf{e}}(k)=\bm{\varsigma}(k)-{{\bm{\varsigma}}_{\bf{L}}}(k)$, and ${{\bm{\varsigma}}_{\bf{R}}}(k)$ denotes the reference image feature without the predictive error. Thus, we can obtain:

$${{\bm{\varsigma}}_{\bf{R}}}(k)-{{\bm{\varsigma}}_{\bf{L}}}(k)={{\bm{\varsigma}}_{\bf{G}}}(k)-\bm{\varsigma}(k)$$

(12)

The object of the visual MPC controller is then transformed into minimizing the tracking error of the prediction model with respect to ${{\bm{\varsigma}}_{\bf{R}}}(k)$. Let ${{\bm{\varsigma}}_{\bf{L}}}(k)={{\bm{\varsigma}}_{\bf{L}}}(k+1)-{{\bm{\varsigma}}_{\bf{L}}}(k)$, Eq. (8) can be rewritten as the following state-space representation:

$$\left\{\begin{matrix}{\bf{x}}(k+1)={\bf{x}}(k)+{\bf{B}}(k){\bf{u}}(k)\\ {\bf{y}}(k)={\bf{x}}(k)\end{matrix}\right.$$

(13)

where the system state ${\bf{x}}(k)={{\bm{\varsigma}}_{\bf{L}}}(k)$, the control variable ${\bf{u}}(k)=\Delta{\bf{q}}(k)$, ${\bf{y}}(k)$ is the output and ${\bf{B}}(k)={\bf{\hat{L_{m}}}}(k){\bf{\hat{J}}}(k)$.

In addition, some constraints should be considered. To ensure that the MicroNeuro remains stable and avoid undesirable contact with the brain ventricles, the camera position should meet certain constraint:

$${\bf{P}}^{\min}\leq{\bf{P}}(k)\leq{\bf{P}}^{\max}$$

(14)

Correspondingly, considering some physical hard constraints on MicroNeuro, such as the restriction of the capability of the motors, actuator constraint is defined as follows:

$${\bf{q}}^{\min}\leq{\bf{q}}(k)\leq{\bf{q}}^{\max}$$

(15)

Moreover, to ensure that the target of concern is always within the field of view and away from areas with large camera distortion, output constrain is described as follows:

$${{\bm{\varsigma}}_{\bf{L}}}^{\min}\leq{{\bm{\varsigma}}_{\bf{L}}}(k)\leq{{\bm{\varsigma}}_{\bf{L}}}^{\max}$$

(16)

At each sample time $k$, the current measured system state is set as the initial state of an optimal control problem (OCP) with constrains, and the current control action is determined by solving the problem in the further $N_{P}$ sampling periods. Only the first optimal input is applied on the system in the optimal input sequence of length $N_{c}$. $N_{p}$ and $N_{c}$ are identified as the prediction horizon and control horizon, respectively. The objective is described as follows:

$$\begin{split}&\min_{{\bf{U}}(k)}{\bf{V}}({\bf{U}}(k))=\sum_{i=k}^{k+N_{p}-1}||{\bf{Y}}(k)-{\bf{S}}(k)||^{2}_{{\bf{Q}}}\\ &=\sum_{i=k}^{k+N_{p}-1}({\bf{y}}(i|k)-{{\bm{\varsigma}}_{\bf{R}}}(i|k))^{\mathrm{T}}{\bf{Q}}({\bf{y}}(i|k)-{{\bm{\varsigma}}_{\bf{R}}}(i|k))\end{split}$$

(17)

subject to Eq. (12), (14), (15) and (16). In Eq. (17), ${\bf{U}}(k)\in\mathbb{R}^{8N_{p}\times 1}$, ${\bf{U}}(k)=({\bf{u}}(k|k)\cdots{\bf{u}}(k+N_{c}-1|k)\cdots{\bf{u}}(k+N_{c}-1|k))^{\mathrm{T}}$ is the control sequence, ${\bf{Y}}(k),{\bf{S}}(k)\in{\mathbb{R}}^{2N_{p}\times 1}$ are output and reference sequence, ${\bf{Q}}\in\mathbb{R}^{2\times 2}$ is the weight matrix. ${\bf{y}}(i|k)$ denotes the predictive value of output at $i$-th sample time. Problem (17) can further come down to solve a quadratic programming (QP) with constrains. Specially, in our implementation, problem (17) is formulated in CasADi [30] and is solved using its built-in optimization solvers.

In this experiment, the region of interest (ROI) was defined as the center of the image ${{\bm{\varsigma}}_{\bf{G}}}(k)=(355,355)^{\mathrm{T}}$. As shown in Fig. 5(b), the MicroNeuro system was commanded to bring the specifically chosen markers to the ROI, which were distributed at $60^{\circ}$ intervals on a printed circle. The experimental analysis involved conducting six trials, and the effectiveness was demonstrated through the measured trajectories of the markers, as depicted in Fig. 5(c). In each instance, the robot successfully returned the marker to the center with average terminal error was 21.8 pixels. The average time required to complete the tracking task across the six experiments was measured to be 11.25 s. This accomplishment highlights the robustness and reliability of the proposed method in achieving fast and precise tracking.

The experiment was designed to evaluate the stability of the proposed system following a target in a dynamic environment. As shown in Fig. 6(a), the robot tracked an AprilTags marker attached to a linear guide, positioned 20mm from the robot’s camera. The guide reciprocated at a speed of $2.5mm/s$ over a $20mm$ stroke. Fig. 7 illustrates that tracking errors decreased significantly once the marker was captured by the camera, with errors reduced to below the MPE within 6 s in Test 1, reaching a lowest error of 2.23 pixels. After the initial stable tracking of the target was accomplished, the standard deviation (SD) of the errors for test 1 and 2 were 20.85 and 21.81 pixels, respectively, which further supports the effectiveness of the system in maintaining precise tracking of the target.

This experiment was designed to evaluate the robot’s ability to follow a set trajectory that guides the marker along a defined path in the captured image. Experiment setup was same as Fig. 5(a). Under the guidance of the controller, the robot automatically completes tracking of multiple key target points on different trajectories to approximately complete the tracking of curves in the image plane. These discrete key target points set on the letters ${CAIR}$. The experimental results in Fig. 8 showed that the controller has good tracking performance for the key points of each trajectory. The root mean square error (RMSE) of the four curves were 11.66, 11.62, 11.30 and 11.95 pixels respectively.

In clinical procedures, the use of endoscopic instruments like biopsy gripper and electrocoagulation, inserted via the working channel, can significantly disrupt the flexible endoscope’s view, leading to loss of lesion visibility or inadequate operating angles. This experiment aims to assess the robustness of the proposed method against external disturbances, ensuring the endoscope stays focused on the ROI. In the 3D printed brain shown in Fig. 6(b), we placed a marker in the pineal gland region to mark the area of interest. Initially, the robot was manually operated to roughly approach the target area through one burr hole, and the visual MPC controller has quickly tracked the target, as shown in Fig. 9. The insertion and operation of biopsy forceps introduced rapid noise to the robot, significantly increasing tracking error. However, the controller adjusted the tool within ten steps, reducing the error to less than 30 pixels. This result demonstrates the controller’s ability to enhance the MicroNeuro robot’s resistance to interference, suggesting its potential application in neurosurgery.