Bimanual crop manipulation for human-inspired robotic harvesting

Sotiris Stavridis^*, Dimitrios Papageorgiou, Leonidas Droukas and Zoe Doulgeri
The research leading to these results has received funding from the European Community’s Framework Programme Horizon 2020 under grant agreement No 871704, project BACCHUS.
Authors are with the Automation & Robotics Lab, School of Electrical & Computer Engineering, Aristotle University of Thessaloniki, Greece. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>^*Corresponding Author

Most existing robotic harvesters utilize a unimanual approach; a single arm grasps the crop and detaches it, either via a  detachment movement, or by cutting its stem with a specially designed gripper/cutter end-effector.
However, such unimanual solutions cannot be applied for sensitive crops and cluttered environments like grapes and a vineyard where obstacles may occlude the stem and leave no space for the cutter's placement.  In such cases, the solution would require a bimanual robot in order to visually unveil the stem and manipulate the grasped crop to create cutting affordances which is similar to the practice used by  humans.
In this work, a dual-arm coordinated motion control methodology for reaching a stem pre-cut state is proposed. The camera equipped arm with the cutter is reaching the stem, unveiling it as much as possible, while the second arm is moving the grasped crop towards the surrounding free-space to facilitate its stem cutting. Lab experimentation on a mock-up vine setup with a plastic grape cluster evaluates the proposed methodology, involving two UR5e robotic arms and a RealSense D415 camera.

§ INTRODUCTION

In the last decade, development of robotic technologies and their application in agricultural tasks are becoming a growing topic of interest [1]. The rising of food supply demands, climate change, land degradation and arable land limitations have all pushed towards agricultural productivity growth becoming an important priority. Incorporating advanced, automated, agricultural technologies is bound to greatly benefit productivity and in turn economic development [2], while at the same time facilitating and improving the difficult working conditions of farmers and agricultural workers [3]. In this context, an increasing amount of research work has been noticed during the last years [4, 5]. Various agricultural tasks are considered as potential applications of robotic solutions developed in literature e.g. yield estimation, phenotyping, pruning, as well as crop harvesting. Earlier solutions in crop harvesting involved the bulk concept (i.e. tree trunk/branch shaking), however the selective concept has been mostly adopted lately, since it entails considerably less danger of harming the crops. In selective harvesting, the robotic system firstly detects a specific harvest target and then harvests it.

A variety of integrated robotic solutions have been proposed towards this end [6], with the majority of them utilizing one robotic manipulator (unimanual setup), responsible for monitoring the target crop via a mounted in-hand camera and harvesting it, i.e. grasping it and detaching it. In [7], [8] a detaching movement performed by the arm's gripper e.g. pulling and twisting the grasped crop in order to cause the stem's breaking and crop's release. However, detaching movements are only applicable on rigid crops (e.g. apples) that can withstand such treatment without being damaged. Alternatively, various approaches utilize the grasped crops stem cutting and mainly involve: (a) a cutting tool embodied in the gripper that cuts the stem after the crop is grasped [9], [10]; (b) a specially designed cutter end-effector, responsible for simultaneously cutting and holding the crop by its stem [11], [12]. In these works, harvesting of eggplants, cucumbers, tomatoes and strawberries is performed. Such unimanual solution may work after certain crop modification e.g. pruning or cultivation/farming methodologies specifically designed to be more "friendly" towards automated/robotic harvesting such as the cucumber high-wire system [10] and the apple V-trellis system [8]. However, in general, such unimanual solutions cannot be applied for sensitive and cluttered environments like a vineyard where poles, wires, flexible branches and leaves and even the crop itself partially occlude the stem, which can be further wrapped partially around a branch leaving no space for the cutter placement. In such cases the solution would require a bimanual robot in order to visually unveil the stem and manipulate the grasped crop to create cutting affordances. Notice that grape harvesting is performed bimanually by humans. In fact, a very common approach that humans apply during grape harvesting, is to grasp with one arm the crop, moving their head and body to be able to see the stem and manipulate the crop, if needed, in order to provide sufficient free room around the stem for the cutting tool held by their other arm (Fig. <ref>).

Bimanual human harvesting

There are few works in literature, proposing a dual-arm/bimanual solution [13, 14]. In [13], a
dual-arm tomato harvesting robot is designed, with two 3 degrees of freedom (DOF) SCARA manipulators (one prismatic and two revolute joints) and two different end-effectors (vacuum suction gripper and saw cutting tool). In [14], the same robotic solution is further developed, introducing autonomous ripe tomato detection in contrast to [13] where human collaboration was required, evaluated in a simple environment with several potted tomatoes placed in rows. However, the limited manipulability of the arms due to few DOF and the very specific gripper's design considering only rigid, spherical crops would render the above solution's application very difficult in more challenging crops/environments, such as grapes/vineyards.

In this work, we consider a bimanual robot with one arm equipped with a crop specific gripper (grasping arm) and the other with a stem cutting tool and a camera mounted on its end-effector (camera arm). Our focus is on the core bimanual task of approaching while unveiling the stem by the camera arm, while manipulating the grasped crop with the other arm for generating cutting affordances. Reaching to grasp the crop and reaching the final cut pose to cut the stem from a precut pose are mainly unimanual operations, which can be addressed by methods already proposed in the literature, like grasp planning [15] and visual-servoing [16].

§ PROBLEM DESCRIPTION

We consider a robot with two robotic arms: (a) the grasping arm of $N_g \in \mathbb{N}$ degrees of freedom (DOF) with $\vct{q}_g \in \mathbb{R}^{N_g}$ its joints vector, responsible for the stable grasping and manipulation of target crop; (b) the camera arm of $N_c \in \mathbb{N}$ DOF and joint vector $\vct{q}_c \in \mathbb{R}^{N_c}$ that monitors the whole scene via a mounted, in-hand, RGB-D camera equipped with a cutting tool. As an initial state for the robot we consider that the grasping arm has achieved grasping of the crop by its gripper and that a point-cloud of the overall scene is provided by the camera arm. We assume that an appropriate vision algorithm/methodology designed specifically for recognizing the stem is available. Notice that this assumption does not exclude the case of a stem being partially occluded or initially not recognized.

Let $\vct{T}_c \in SE(3)$ be the homogeneous transformation expressing the generalized pose of the camera's frame $\{C\}$ (which is in our case placed at the camera arm's tip) with respect to the world inertial frame $\{0\}$; it involves its position $\vct{p}_c(\vct{q}_c) \in \mathbb{R}^3$ and its orientation $\vct{R}_c(\vct{q}_c) \triangleq [\vct{x}_c \; \vct{y}_c \; \vct{z}_c ] \in SO(3)$, that are both known via the camera arm's kinematics. Additionally, let $\vct{T}_g \in SE(3)$ be the homogeneous transformation expressing the grasping arm's tip/grasped crop frame $\{G\}$ with respect to $\{0\}$, involving position $\vct{p}_g(\vct{q}_g) \in \mathbb{R}^3$ and orientation $\vct{R}_g(\vct{q}_g) \triangleq [\vct{x}_g \; \vct{y}_g \; \vct{z}_g ] \in SO(3)$, known via the grasping arm's kinematics.

The core bimanual task of approaching and unveiling the stem with the camera arm and manipulating the grasped crop with the grasping arm, which is addressed in this work, is a bimanual asymmetric task where each arm's controller has a different objective but their motion needs appropriate planning and coordination. In particular, the camera arm's controller should operate in a clutter dynamical scene in which changes are not only induced by dynamic changes in the environment (e.g wind gusts can move the surrounding leaves) but also by the grasping arm's manipulation of the crop. In the following two subsections each arm's control objectives are presented.

§.§ Camera arm control objectives

The camera arm is responsible for reaching a predefined region of interest (ROI) surrounding the crop and its stem and adopt a pose in this region that allows the stem to be unveiled as much as possible thus getting sufficiently close to it into a pre-cutting position. We can therefore distinguish two objectives that should be simultaneously satisfied by this controller, reaching and centering, and unveiling.

We model the ROI by a sphere and we call the associated stem, object of interest (OOI); let $\vct{p}_r \in \mathbb{R}^3$ be the center of ROI in the inertia frame with $r \in \mathbb{R}_{>0}$ its radius. The reaching and centering objective aims at moving the camera arm's end-effector to any of the ROI points, while aligning the camera's $Z$ axis with the vector pointing to the ROI center. The ROI center can be placed on a point on the stem if it is initially partially recognized or on the grape otherwise. As the camera moves closer and visual feedback provides better estimates of the point-cloud of the grape and its stem it is updated. To prevent discontinuities in the arm's operation the center of ROI displacement is smoothed by a filter. The two objectives can be formulated as follows:

*  Reaching and centering:

for reaching the following should hold true: $\|\vct{p}_r - \vct{p}_c\| \rightarrow\overline{r} \leq r$ for $t\rightarrow\infty$ and some scalar positive $\overline{r}$; regarding centering, the angle between $\vct{p}_r - \vct{p}_c$ and the camera's view direction $\vct{z}_c$ defined as $\theta(\vct{p}_c, \vct{p}_r) \triangleq \text{cos}^{-1}\left(\frac{\vct{z}_c^\intercal(\vct{p}_r - \vct{p}_c)}{\|\vct{p}_r - \vct{p}_c\|}\right)$ should satisfy $\theta(\vct{p}_c, \vct{p}_r)\rightarrow 0$ for $t\rightarrow\infty$.

*  Unveiling: maximize the visible part of the OOI, i.e. maximize the perceived number of points of the point-cloud that belong to the stem.

Notice that $\theta \rightarrow 0$ involved in the first objective maintains the ROI center at the center of the camera's field of view, thus securing the best ROI and OOI viewpoint. Moreover, notice that $\theta$ cannot be more than $\nicefrac{\pi}{2}$, since the ROI center $\vct{p}_r$ is within the camera's field of view, i.e. a pyramid.

Clearly, the reaching and centering objective is achieved by reaching any position in ROI which is associated with two specific orientation dof for centering.  The unveiling objective (<ref>) selects the position in ROI that visually unveils the OOI as much as possible. Hence, all the above objectives correspond to 5 dof with the rotation around the z camera axis denoting a task redundancy of 1 dof. The latter can be specified for properly rotating the cutter wrt the stem.

§.§ Grasping Arm control objective

Regarding the grasping arm's task, the goal is to maximize the free room around the stem  to enable the placement of the cutter. This is achieved by manipulating the grasped grape so that the stem is stretched and moved towards a point in the free space that is most distant from any surrounding obstacle identified from the provided scene's point-cloud. Let this point's position be  $\vct{p}_{gd}$. We assume that  $\vct{f} \in \mathbb{R}^3$, the actual applied force vector at grasping arm's tip, is measurable. Notice that the grape is an object  hinged to the branch by its stem. We can estimate the stem basis from the point-cloud and let $\vct{n}_c \in \mathbb{R}^3$ be the direction from the stem basis to $\vct{p}_{gd}$. We can then formulate the control objective as a force position control objective along $\vct{n}_c$ and its orthogonal subspace respectively.  Hence if  $f_d \in \mathbb{R}$ is an appropriately selected force magnitude for stretching the stem sufficiently without damaging/tearing it, the grasping arm's control objectives can be mathematically expressed as follows: $\vct{n}_c \vct{n}_c^\intercal \vct{f} \rightarrow \vct{n}_c f_d$ and $\vct{p}_{g} \rightarrow (\vct{I}_3-\vct{n}_c \vct{n}_c^\intercal)\vct{p}_{gd}$.

Notice that we have decomposed the space of the grasping arm's end-effector position  based on the identified $\vct{p}_{gd}$ and thus both  the position and force control action space can create motion until the final target is reached. Alternatively the space should be decomposed appropriately given the current position of the grape and its stem. As compared to the alternative decomposition solution, the adopted approach is computationally lighter and shown to be effective as we will demonstrate in the experimental part.

§ PROPOSED CONTROL METHODOLOGY

We consider a velocity controlled bimanual robot. Hence our control objective is to design a reference velocity control signal $\vct{V}_r = [\vct{V}_c^\intercal \ \ \vct{V}_g^\intercal]^\intercal \in \mathbb{R}^{12}$ involving the arm end-effector velocities $\vct{V}_c, \ \vct{V}_g \in \mathbb{R}^6$ for the cutting and grasping arm respectively, which can then be mapped to the joint space via the robot's extended Jacobian matrix $\vct{J}(\vct{q}) = diag(\vct{J}_c(\vct{q}_c) , \vct{J}_g(\vct{q}_g)) \in \mathbb{R}^{12 \times (N_c + N_g)}$ involving the camera arm Jacobian $\vct{J}_c(\vct{q}_c) \in \mathbb{R}^{N_c \times 6}$ and the grasping arm Jacobian $\vct{J}_g(\vct{q}_g) \in \mathbb{R}^{N_g \times 6}$  as follows:

\begin{equation} \label{eq:controlled_system}
\dot{\vct{q}} = \vct{J}^\dagger(\vct{q}) \vct{V}_r .
\end{equation}

where  $\vct{q} = [\vct{q}_c^\intercal \ \ \vct{q}_g^\intercal]^\intercal \in \mathbb{R}^{N_c + N_g}$ denotes the robot's joint position, $\vct{J}^\dagger(\vct{q})$ is a pseudo inverse of the robot's extended Jacobian matrix.

Designing $\vct{V}_c$ and $\vct{V}_g$ is heavily based upon the overall scene's point-cloud set and its processing in order to estimate the required quantities involved in each control objective as well as to identify the point-cloud of the stem and the surrounding obstacles.

§.§ Scene Point-cloud Process and   calculations of control related critical point positions

Scene's point-cloud $\mathcal{W}$. Obstacles subset $\mathcal{O}$: blue. Stem subset $\mathcal{S}$: cyan.

Stem's point-cloud clustering and line segment fitting: green and red. Stem's base: blue x.

Obstacles projection upon sphere subset $\mathcal{O}_{pr}$: purple.

Sphere surface sampling ($500$ samples). Near-obstacles sampled points: black. $\mathcal{O}$ and $\mathcal{S}$ omitted for clarity.

Free-space sphere points: green. Near-obstacles sampled points: black. $\mathcal{O}_{pr}$ omitted for clarity.

Free-space hemisphere subset $\mathcal{F}$: green. Obstacles $\mathcal{O}$: blue. Plane fitted into $\mathcal{O}$: black lines. Near-obstacles sampled points omitted for clarity.

Point-cloud processing example from lab experimentation with a Real-Sense camera and a mock-up crop. Figures shown from camera's viewpoint with $Z$ axis corresponding to the camera's view direction.

Let $\mathcal{W}: \{ \vct{p}_{w_i} \in \mathbb{R}^3 \ , \ i=1,...,n_W \}$ be the set of $n_W \in \mathbb{N}$ points of the point-cloud captured by the camera  with $\vct{p}_{w_i}$ being their position vector  expressed in the world frame (Fig. <ref>). Let subset $\mathcal{S} \subseteq \mathcal{W}: \{ \vct{p}_{s_j} \in \mathbb{R}^3 \ , \ j=1,...,n_S \}$ denote the detected stem's point-cloud of $n_S \in \mathbb{N}$ points of position $\vct{p}_{s_j}$ (Fig. <ref>, cyan points). To identify the stem's base $\vct{p}_{sb} \in \mathbb{R}^3$, we model the stem by connected linear segments to account for their bending capability. To this end we initially partition $\mathcal{S}$ into separate, consecutive point-cloud clusters, utilizing any of the existing clustering methodologies, e.g. the k-means clustering algorithm. Since crop stems are generally short (e.g. $5-7$ cm), two clusters are considered to be enough to adequately model the stem's flexibility. We then identify which cluster is further and closer to the grasping arm, comparing the distance between the gripper and each cluster mean value and denote them as "top" $\mathcal{T}: \{ \vct{p}_{cl_a} \in \mathbb{R}^3 \ , \ a=1,...,n_T \}$ and "bottom" $\mathcal{B}: \{ \vct{p}_{cl_b} \in \mathbb{R}^3 \ , \ b=1,...,n_B \}$ cluster respectively where $\vct{p}_{cl_a}, \ \vct{p}_{cl_b}$ and $n_T,n_B$ are their respective points' position vectors and number, with $\mathcal{T},\mathcal{B} \subseteq \mathcal{S}$ and $\mathcal{T} \neq \mathcal{B}$. Utilizing the Principal Component Analysis (PCA), a line segment is fitted into each cluster. Let  $\vct{n}_{st}, \ \vct{n}_{sb} \in \mathbb{R}^3$ be their direction vectors found by the corresponding major PCA eigenvectors. We can also find the line segments'  length $l_t, \ l_b \in \mathbb{R}_{>0}$ for the top and bottom line respectively, and deduce an estimate of the whole stem's length $l = l_t + l_b \in \mathbb{R}_{>0}$.

Point-cloud Subsets     Involved Points
$\text{whole scene}: \mathcal{W}$
$\vct{p}_{w_i}  \in \mathbb{R}^3 \ , \ i=1,...,n_W$

$\text{stem}:        \mathcal{S} \ , \ \mathcal{S} \subseteq \mathcal{W}$
$\vct{p}_{s_j}  \in \mathbb{R}^3 \ , \ j=1,...,n_S$

$\text{stem top}:    \mathcal{T} \ , \ \mathcal{T} \subseteq \mathcal{S}$
$\vct{p}_{cl_a} \in \mathbb{R}^3 \ , \ a=1,...,n_T$

$\text{stem bottom}: \mathcal{B} \neq \mathcal{T} \ , \ \mathcal{B} \subseteq \mathcal{S}$
$\vct{p}_{cl_b} \in \mathbb{R}^3 \ , \ b=1,...,n_B$

$\text{obstacles}:   \makecell[l]{\mathcal{O} \triangleq (\mathcal{W} - \mathcal{S}) \cap \\ S_{ph}(\vct{p}_{sb},r_o)}$
$\vct{p}_{o_k}  \in \mathbb{R}^3 \ , \ k=1,...,n_O$

$\text{proj. obstacles}: \makecell[l]{\mathcal{O}_{pr} \triangleq \textit{proj}(\mathcal{O}) \ \\ \text{upon} \ \text{S}(\vct{p}_{sb},l)}$
$\textit{proj}_{\text{S}}(\vct{p}_{o_k})  \in \mathbb{R}^3 \ , \ k=1,...,n_O$

$\text{free-space}:  \mathcal{F}$
$\vct{p}_{f_m}  \in \mathbb{R}^3 \ , \ m=1,...,n_F$

Point-cloud Classification.

The stem's base $\vct{p}_{sb}$ is then calculated as follows (Fig. <ref>):

\begin{align} \label{stem_basend}
\vct{p}_{sb} &= \text{mean}(\mathcal{T}) - (l_t/2)\vct{n}_{st}
\end{align}

where $\text{mean}(X)$ is the mean value of $X$ point-cloud subset. For simplicity, we consider that target crop's stem is attached to a branch with limited mobility. Hence, the stem's base $\vct{p}_{sb}$ is assumed to be static and is identified only once at the start of the reaching/unveiling task.

All scene points that do not belong to $\mathcal{S}$ can be considered as potential obstacles but we limit our search for obstacles in an area closely around the stem, i.e. a sphere $\text{S}(\vct{p}_{sb},r_o)$ with the stem's base $\vct{p}_{sb}$ as its center and a radius of $r_o \in \mathbb{R}$, selected freely to be adequately close to the stem. Thus, the obstacle point-cloud subset consists of $n_O \in \mathbb{N}$ points, that are not stem points and reside inside the above considered sphere, with vector $\vct{p}_{o_k}$ describing their position: $\mathcal{O}: \{ \vct{p}_{o_k} \in \mathbb{R}^3 \ , \ k=1,...,n_O \}$, with $\mathcal{O} \triangleq (\mathcal{W} - \mathcal{S}) \ \cap \ \text{S}(\vct{p}_{sb},r_o)$ (Fig. <ref>, blue points).

Let subset $\mathcal{F}: \{ \vct{p}_{f_m} \in \mathbb{R}^3 \ , \ m=1,...,n_F \}$ denote the free-space surrounding the stem, containing $n_F \in \mathbb{N}$ points with position vector $\vct{p}_{f_m}$. In order to calculate $\mathcal{F}$, we define a sphere $\text{S}(\vct{p}_{sb},l)$ with the stem's base $\vct{p}_{sb}$ as its center and a radius equal to the calculated stem's length $l$. Subset $\mathcal{O}$ points are then projected upon this sphere's surface with $\mathcal{O}_{pr}: \{ \textit{proj}_{\text{S}}(\vct{p}_{o_k}) \in \mathbb{R}^3 \ , \ k=1,...,n_O \}$ denoting the set of these projected points with $\textit{proj}_{\text{S}}(\vct{p}_{o_k})$ their position vector (Fig. <ref> - <ref>, purple points) given by: $\textit{proj}_{\text{S}}(\vct{p}_{o_k}) = \frac{\vct{p}_{o_k} - \vct{p}_{sb}}{\|\vct{p}_{o_k} - \vct{p}_{sb}\|}l + \vct{p}_{sb}$. Sphere $\text{S}(\vct{p}_{sb},l)$ is uniformly sampled utilizing the Fibonacci lattice methodology, which is a fast and near optimal way of sampling a sphere. Furthermore, in [17] a method to find the closest Fibonacci point for a given point upon the sphere is introduced, allowing the match of the sphere projected obstacle points $\mathcal{O}_{pr}$ to corresponding Fibonacci points without resorting to an exhaustive search that would significantly increase the calculation cost. Hence assuming adequate sampling, the sampled sphere points nearest to $\mathcal{O}_{pr}$ points are found and then removed from the sphere's surface, leaving the remaining sphere surface points which constitute the free-space surrounding the crop's stem (Fig. <ref>).
Since the camera cannot see behind obstacles (e.g. leaves, target crop), we apply PCA on subset $\mathcal{O}$ to calculate the corresponding two major eigenvectors i.e. the two direction vectors where most $\mathcal{O}$ data reside; this two vectors define a representative plane fitted into $\mathcal{O}$ (Fig. <ref>, black lines). Any free-space point upon the sphere's surface which is behind this plane, is considered invalid. Hence, the previously defined sphere becomes a hemisphere, ensuring that the determined free-space $\mathcal{F}$ is deduced so that it always resides within the camera's clearly visible field of view (Fig. <ref>).
From $\mathcal{F}$, free-space point $\vct{p}_{gd}$, which is tied to the grasping arm's control objective (Subsection <ref>), is calculated by solving the following optimization problem:

\begin{equation} \label{eq:attr_point}
\vct{p}_{gd} = \operatorname*{argmax}_{\vct{p}_{f_m}} \Big\{ \min_{\forall k} \big\{ \|\vct{p}_{f_m} - \textit{proj}_{\text{S}}(\vct{p}_{o_k})\| \big\} \Big\},
\end{equation}

utilizing the projected obstacle points $\textit{proj}_{\text{S}}(\vct{p}_{o_k})$ upon the sphere $\text{S}(\vct{p}_{sb},l)$. Given a relatively low number of free-space points by appropriate selection of the sphere's sampling rate, (<ref>) can be solved with brute force in a short amount of time. Table <ref> summarizes the classified/calculated point-cloud subsets.

§.§ Camera arm Control - Reaching & Unveiling Tasks

The camera arm end-effector's reference velocity is calculated by the superposition of two reference velocity control terms, the ROI reaching with centering control signal $\vct{V}_{cr} \in \mathbb{R}^6$ and the  OOI unveiling control signal, $\vct{V}_{cu} \in \mathbb{R}^6$:

\begin{equation} \label{eq:superimposition}
\vct{V}_c = \vct{V}_{cr} + \vct{V}_{cu},
\end{equation}

When there are no stem points detected ($n_S = 0$) $\vct{V}_{cu} = \vct{0}_{6 \times 1}$  by design and  $\vct{V}_{cr}$ is acting alone. Such a case may occur during the start of the motion, when the camera is far from the ROI and the OOI starts being visible only when it approaches or enters the ROI. Figure <ref> depicts a typical  initial and final state of the camera's motion under the superimposed velocities.

§.§.§ Reaching reference velocity

Mathematically, reaching is fulfilled when the camera position $\vct{p}_c$ converges to the manifold $\Omega \triangleq \{ \vct{p}_c \in \mathbb{R}^3: \text{f}(\vct{p}_r - \vct{p}_c) \leq 0 \}$, with $\text{f}(\vct{x}) \triangleq \vct{x}^\intercal \vct{x} - r^2$ for any $\vct x \in \mathbb{R}^3$. Centering the camera's field of view with respect to ROI center regards the orientation of the camera and is achieved when the camera's $Z$ axis points to the center of the ROI (Fig. <ref>). To achieve the reaching with centering objective we propose the following control signal, which draws its inspiration from the region reaching approach proposed in [18] for the translational part:

\begin{equation} \label{eq:control_signal_rr_world}
\vct{V}_{cr} \triangleq
\begin{bmatrix}
k_{cp} \text{max}(0,\; \text{f}(\vct{e}_{pc}))\vct{e}_{pc} \\
k_{co} \theta \vct{k}
\end{bmatrix},
\end{equation}

where $k_{cp}, k_{co} \in\mathbb{R}_{>0}$ are constant positive gains for translation and orientation respectively and $\vct{e}_{pc} = \vct{p}_r - \vct{p}_c$; $\theta\in[0,\nicefrac{\pi}{2}]$ and $\vct k$ are the angle and axis of the minimum rotation between $\vct{e}_{pc}$ and $\vct{z}_c$, which can be calculated by the following expressions:

\begin{equation} \label{eq:ktheta_world}
\theta = \text{cos}^{-1} \left( \frac{\vct{z}_c^\intercal \vct{e}_{pc}}{\|\vct{e}_{pc}\|}\right), \;
\vct{k} = \frac{\vct{S}( \vct{z}_c ) \vct{e}_{pc}}{\|\vct{S}( \vct{z}_c )\vct{e}_{pc}\|}
\end{equation}

with $\vct{S}(\vct{z})$ the skew symmetric matrix of the corresponding vector $\vct{z}$. Notice that $\vct{V}_{cr}$ is continuous and its value is zero if and only if $\text{f}(\vct{e}_{pc}) \leq 0$ and $\theta = 0$. When this signal is acting alone it is easy to prove asymptotic convergence to the manifold $\text{f}(\vct{e}_{pc}) \leq 0$ and $\theta = 0$.

§.§.§ Unveiling Reference velocity

Regarding unveiling, we consider obstacles not just as points but as spheres with radius $d_o \in \mathbb{R}$ selected such that the empty space between neighboring points of the obstacle point-cloud subset $\mathcal{O}$ is covered (Fig. <ref>). Notice the ray from the camera  to  a visible point of the OOI, $\vct{p}_{s_j}$, shown in Fig. <ref>. Let $\hat{\vct{p}}_{j,k}$ be a point in this ray that has a minimum distance from an obstacle point $\vct{p}_{o_k}$. Let $\hat{r}_{j,k} \in \mathbb{R}_{\geq 0}$ be the distance from the obstacle surface i.e.  $\hat{r}_{j,k} \triangleq \|\hat{\vct{p}}_{j,k} - \vct{p}_{o_k}\| - d_o$. The basic concept of the unveiling reference velocity is targeting at rotating the camera around $\vct{p}_{s_j}$ in order to guarantee that $\hat{r}_{j,k}$ will increase and $\hat{r}_{j,k} > d_o$ so that the unveiled region around $\vct{p}_{s_j}$ will also increase.

Camera arm control methodology. (a) Initial state. (b) Desired state.

Let $\vct{p}_{s_j} \ , \ j=1,...,n_S$ be the currently visible OOI points and $\vct{p}_{o_k} \ , \ k=1,...,n_O$ obstacles.
The nearest point of the $j$-th ray from the $k$-th obstacle is calculated as follows:

\begin{equation}\label{eq:nearest}
\hat{\vct{p}}_{j,k} \triangleq \begin{cases}
\vct{p}_c,\; &\text{if} \; \Big( \frac{\vct{\bar{p}}_{cs_j}^\intercal}{\|\vct{\bar{p}}_{cs_j}\|} \Big) \vct{\bar{p}}_{co_k} \leq 0 \\

\vct{p}_c + \frac{\vct{\bar{p}}_{cs_j}\vct{\bar{p}}_{cs_j}^\intercal}{\vct{\bar{p}}_{cs_j}^\intercal\vct{\bar{p}}_{cs_j}} \vct{\bar{p}}_{co_k}, \; &\text{if} \; \Big( \frac{\vct{\bar{p}}_{cs_j}^\intercal}{\|\vct{\bar{p}}_{cs_j}\|} \Big) \vct{\bar{p}}_{co_k} \in(0, \|\vct{\bar{p}}_{cs_j}\|) \\

\vct{p}_{s_j}, \; &\text{if} \; \Big( \frac{\vct{\bar{p}}_{cs_j}^\intercal}{\|\vct{\bar{p}}_{cs_j}\|} \Big) \vct{\bar{p}}_{co_k} \geq \|\vct{\bar{p}}_{cs_j}\| \\
\end{cases}
\end{equation}

with $\vct{\bar{p}}_{cs_j} = \vct{p}_{s_j} - \vct{p}_{c}$ and $\vct{\bar{p}}_{co_k} = \vct{p}_{o_k} - \vct{p}_{c}$.

Drawing our inspiration from our previous work [19], we propose the utilization of a barrier artificial potential field around each obstacle, which induces a virtual repulsive velocity $\vct{u}_{j,k}$, acting at $\hat{\vct{p}}_{j,k}$.
The proposed barrier artificial potential function
is designed to induce a repulsive velocity only within a predefined distance $d_a \in \mathbb{R}_{>0}$ from the obstacle surface and is defined as follows:

\begin{equation}\label{eq:art_potential}
\begin{cases}
\frac{1}{2}\ln^2\left(\frac{d_a^2}{d_a^2-(d_a-\hat{r}_{j,k})^2}\right), \; &\text{if} \; \hat{r}_{j,k}< d_a \\
0, \; & \text{otherwise}
\end{cases},
\end{equation}

Notice that the value of the artificial potential tends to infinity, when the surface of the obstacle is reached from outside in order to guarantee obstacle/occlusion avoidance.

The virtual commanded repulsive velocity is then given by:

\begin{equation}\label{eq:u_ij}
\begin{split}
&\vct{u}_{j,k} \triangleq - \frac{\partial V(\hat{\vct{p}}_{j,k})}{\partial\hat{\vct{p}}_{j,k}} = \\
\frac{2}{d_a^2-(d_a-\hat{r}_{j,k})^2}\ln\left(\frac{d_a^2}{d_a^2-(d_a-\hat{r}_{j,k})^2}\right)\vct{e}_{j,k}, \;  &\text{if} \; \hat{r}_{j,k}< d_a \\
\vct{0}_{3 \times 1}, \; &\text{otherwise}
\end{cases}
\end{split}
\end{equation}

with $\vct{e}_{j,k} \triangleq (d_a-\hat{r}_{j,k})\frac{\hat{\vct{p}}_{j,k} - \vct{p}_{o_k}}{\|\hat{\vct{p}}_{j,k} - \vct{p}_{o_k}\|} \in \mathbb{R}^3$. The virtual repulsive velocity (<ref>) possesses the following properties:

* $\|\vct{u}_{j,k}\| \neq 0$, if and only if $0 < \hat{r}_{j,k} < d_a$, and $\vct{u}_{j,k} \neq \vct{0}_{3 \times 1}$ if and only if $\hat{r}_{j,k} \geq d_a$, which means that the signal will not be affected by rays that are not within the range of influence of the $k$-th obstacle point, defined by $d_a$.

* $\|\vct{u}_{j,k}\| \rightarrow \infty$, when $\hat{r}_{j,k} \rightarrow 0$, i.e. when the $j$-th ray approaches the surface of the $k$-th obstacle. Notice that $\|\hat{\vct{p}}_{j,k} - \vct{p}_{o_k}\|$ cannot be less than $d_o$, as it reflects the accuracy of the RGB-D camera, by definition.

* $\vct{u}_{j,k}$ is continuous with respect to $\hat{\vct{p}}_{j,k}$.

* $\vct{u}_{j,k}$ is, in general, not continuous in time, as $\vct{p}_{s_j}, \vct{p}_{o_k}$ and even $n_S, \ n_O$ depend on the point-cloud perceived by the RGB-D camera during its motion; as more points of the OOI, previously occluded, enter into the field of view, it is possible that some of them induce a non-zero control signal. This discontinuity can be remedied by a first order low-pass filter.

To synthesize the total proposed control signal $\vct{V}_{cu}$, we calculate the angular velocity $\vect{\omega}_{j,k}$, which is induced by the virtual repulsive velocity $\vct{u}_{j,k}$ around an axis passing from $\vct{p}_{s_j}$ and defined by the cross product of the directions of $\vct{u}_{j,k}$ and $\vct{p}_{s_j}$. This angular velocity is given by:

\begin{equation}\label{eq:omega_pivot}
\vect{\omega}_{j,k}\triangleq
\begin{cases}
\frac{\vct{S}(\hat{\vct{p}}_{j,k} - \vct{p}_{s_j})}{\|\hat{\vct{p}}_{j,k}-\vct{p}_{s_j}\|^2} \vct{u}_{j,k}, \; &\text{if} \; \Big( \frac{\bar{\vct{p}}_{cs_j}^\intercal}{\|\bar{\vct{p}}_{cs_j}\|} \Big) \bar{\vct{p}}_{co_k} \in(0, \|\bar{\vct{p}}_{cs_j}\|) \\
\vct{0}_{3 \times 1}, \; &\text{otherwise}
\end{cases}.
\end{equation}

By summing the $\vect{\omega}_{j,k}$-s acting on the $j$-th pivot point, we get the total angular velocity for the $j$-th OOI point, which is given by:

\begin{equation}\label{eq:omega_i}
\vect{\omega}_{s_j} \triangleq \sum_{k=1}^{n_O} \vect{\omega}_{j,k} \in \mathbb{R}^3.
\end{equation}

Lastly, to synthesize the total control signal $\vct{V}_{cu}$, the $\vect{\omega}_{s_j}$ are superimposed after calculating the corresponding linear velocity at the end-effector. This superimposition is given by:

\begin{equation}\label{eq:v_c2}
\vct{V}_{cu}
\triangleq k_c \sum_{j=1}^{n_S}
\begin{bmatrix}
\vct{S}(\bar{\vct{p}}_{cs_j}) \\
\vct{I}_3
\end{bmatrix}
\vect{\omega}_{s_j} \in \mathbb{R}^6,
\end{equation}

where $k_c \in \mathbb{R}_{>0}$ is a positive, tunable gain.

Taking the time derivative of the artificial potential (<ref>), we find that: $\dot{V} = - \frac{\|\bar{\vct{p}}_{cs_j}\|-\|\hat{\vct{p}}_{j,k}-\vct{p}_c\|}{\|\vct{p}_{s_j}-\hat{\vct{p}}_{j,k}\|}\vct{u}_{j,k}^\intercal\vct{u}_{j,k}$, which is less or equal than $0$, given that $\|\hat{\vct{p}}_{j,k}-\vct{p}_c\|\leq\|\bar{\vct{p}}_{cs_j}\|$, which is true by construction. This means that $V$ is bounded and that $\hat{r}_{j,k}$ is increasing within the area of influence, due to the fact that $V(\hat{r}_{j,k})$ is a decreasing function of $\hat{r}_{j,k}$. Therefore, $\vct{p}_{s_j}$ will not be occluded by the obstacle centered at $\vct{p}_{o,k}$. Furthermore, given that the linear velocity  of $\vct{V}_{cu}$ is orthogonal to $\bar{\vct{p}}_{cs_j}$, which implies that $\|\bar{\vct{p}}_{cs_j}\|$ remains constant, the maximum radius of visibility around $\vct{p}_{s_j}$ is increased with the proposed control signal. As a result, the progressive unveiling of more points of the OOI occurs in a chain-reaction manner, i.e. by unveiling progressively more and more OOI points.

§.§ Grasping Arm Control - Cutting Affordance Task

To realise the grasping arm's control objective (Subsection <ref>), reference velocity $\vct{V}_g$ is designed as a force/position controller, presented below:

\begin{align}
\vct{v}_p &= - k_{pg} ( \vct{I}_{3 \times 3} - \vct{n}_c \vct{n}_c^\intercal ) \vct{e}_{pg} \ \in \mathbb{R}^3 \label{eq:gr_pos_contr} \\
\vct{v}_f &= - \vct{n}_c \bigg( k_{fp} \vct{e}_{fg} + k_{fi} \int_{0}^{t} \vct{e}_{fg} \bigg)\ \in \mathbb{R}^3 \label{eq:gr_forc_contr},
\end{align}

where $\vct{e}_{pg} = \vct{p}_{g} - \vct{p}_{gd} \in \mathbb{R}^3$, $\vct{e}_{fg} = \vct{n}_c^\intercal f - f_d \in \mathbb{R}$ are the position and force errors, $k_{pg}, k_{fp}, k_{fi} \in \mathbb{R}$ are positive, control gains and

\begin{equation} \label{eq:contr_dir}
\vct{n}_c = \frac{ \vct{p}_{gd} - \vct{p}_{sb} }{ \| \vct{p}_{gd} - \vct{p}_{sb} \| }.
\end{equation}

The force control signal (<ref>) is applied in the direction of $\vct{n}_c$ (<ref>) to maximize the stem's stretching at the end of its motion (Fig. <ref>)  and the position control signal (<ref>) along the space orthogonal to $\vct{n}_c$.  Both signals contribute to the motion of
the grasped crop towards its final position which is determined by the  projection of the space target $\vct{p}_{gd}$ on the space orthogonal to $\vct{n}_c$, while the position on the $\vct{n}_c$ direction is indirectly determined by the achievement of the stretching force.

Regarding the orientation of the grasped crop, let the grasped crop's direction coincide with one of the axes of the grasping arm end-effector's frame $\{G\}$, which is assumed to be the y-axis $\vct{y}_{g} \in \mathbb{R}^3$ (Fig. <ref>). The following proposed signal:

\begin{equation} \label{eq:rot_contr}
\vct{v}_{g\omega} = - k_{og} ( \vct{n}_c \times \vct{y}_{g} ) \in \mathbb{R}^3
\end{equation}

with $k_{og}$ a positive gain, aligns the arm's end-effector/grasped crop's direction $\vct{y}_{g}$ with the stretched stem's $\vct{n}_c$, facilitating further the harvesting process by avoiding any stem obstructions, as well as introducing a more human-like manipulation of the target crop (see Fig. <ref>, where in the final configuration just before the stem's cutting, the human hand along with the grasped crop are aligned with the crop's stem).

Combining (<ref>), (<ref>) and (<ref>), the reference velocity $\vct{V}_g$ regarding the grasping arm's end-effector is:

\begin{equation} \label{eq:superimp_g}
\vct{V}_g = \begin{bmatrix} \vct{v}_{gt} \\ \vct{v}_{g\omega} \end{bmatrix} = \begin{bmatrix} \vct{v}_p + \vct{v}_f \\ \vct{v}_{g\omega} \end{bmatrix}
\end{equation}

with $\vct{v}_{gt} = \vct{v}_p + \vct{v}_f \in \mathbb{R}^3$ the grasping arm's translational velocity.

Grasping arm force/position and orientation control. Force control (<ref>) applied at $\vct{n}_c$. Position control (<ref>) applied at subspace $\vct{I}_{3 \times 3} - \vct{n}_c \vct{n}_c^\intercal$. Orientation control aligns $\vct{n}_c$ with $\vct{y}_g$.

§.§ Bimanual robot motion scheduling and coordination

As the camera arm is initially positioned away from the region of reaching we initiate the reaching/unveiling camera arm motion until some preset thresholds before both arms are provided by their respective reference velocities $\vct{V}_g$ (<ref>) and $\vct{V}_c$ (<ref>) to complete the task. Thresholds regard the arm's velocity and the camera's distance from the surface of the region's sphere with radius $r$. This ensures that a locally optimal view of OOI (i.e. stem) is achieved, and a more accurate ROI center $\vct{p}_{sb}$  and free space target $\vct{p}_{gd}$ can be calculated by the current point-cloud  to enable a successful  fulfillment of the bimanual task.  An enhancement of the bimanual motion is also proposed by superimposing the grasping arm's translational velocity to the camera arm reference velocity to assist in avoiding stem occlusions by the gripper motion.

In particular, the translational part of the grasping arm's velocity $\vct{v}_{gt}$ is superimposed at the camera end-effector's velocity $\vct{V}_c$, while an angular velocity component is additionally applied to avoid $\vct{v}_{gt}$ affecting the stem centering. Therefore, the updated camera arm's velocity is given by:

\begin{equation}\label{eq:final_vc}
\vct{V}_{c}' = \vct{V}_{c} + \begin{bmatrix} \vct{v}_{gt} \\  \vct{S}(\vct{p}_c - \vct{p}_{sb})\vct{v}_{gt} \end{bmatrix}
\end{equation}

Note that the angular part of $\vct{V}_g$ is not applied at the camera arm as it might induce velocities that do not correspond to the stem's velocity.

§ EXPERIMENTAL EVALUATION

For the experimental evaluation of the proposed method, a lab setup was built with a mock-up vine illustrated in Fig. <ref>, involving plastic vine branches and leaves as the surrounding obstacles and a plastic grape cluster as the target crop, with a stem of actual length $l_{act} = 0.07$ m. Identifying the point-cloud belonging to the stem is beyond the scope of this work. Thus for the experiments and in order to have this point-cloud available as the proposed method assumes we have used a red colored stem. Two UR5e robotic arms with a control cycle of $2$ ms were utilized. A 3d printed end-effector was used in the grasping arm to hold stably the mock up grape (Fig. <ref>-(c)). On the camera arm (Fig. <ref>-(b)), a Realsense D415 was mounted to capture the RGB image and the scene's point-cloud, with a frame rate of $30$ FPS and a resolution of $848 \times 480$ pixels and a 3d printed end-effector was attached at its tip to resemble a real cutting tool (Fig. <ref>-(c)). To reduce computational complexity, the point-cloud was down-sampled by a factor of 3.

(a) Lab setup. (b) Camera arm. (c) Grasping arm initial state.

The ROI is defined as a sphere with radius $r = 0.35$ m centered at the stem's base $\vct{p}_{sb}$, identified as in Subsection (<ref>). The stem is reliably detectable only within a range of $0.7$ m from the camera. Any other points are obstacles, for which we use  spheres with radius $d_o = 0.001$ m. The control parameters were selected to be $k_{cp} = k_{co} = 1$ for (<ref>), $d_a = 0.01$ m, $k_c = 0.00025$ for (<ref>). Regarding the thresholds that mark the beginning of the bimanual motion, we use the camera's distance from the ROI sphere a to be  $1.05 r$ and tip translational and angular velocities to be less than $0.01$ m/s and $0.025$ rad/s respectively.

With target $\vct{p}_{gd}$ calculated from (<ref>) and considering as obstacles all non-stem points residing in sphere $\text{S}(\vct{p}_{sb},r_o)$ with $r_o$ equal to the estimated  stem length $l = 0.0631$ m, the grasping arm's control parameters in (<ref>) are chosen to be $k_{pg} = 0.15$, $k_{og} = 0.3$, $k_{fp} = 0.001$, $k_{fi} = 0.0002$, with a desired force magnitude of $f_{d} = 3$ N. Notice that obstacle points are limited within $\text{S}(\vct{p}_{sb},r_o)$ only during the grasping arm's control, since in this case we are mostly interested in obstacles close enough to the stem, surrounding it; in the camera arm's case, we consider all non-stem points as obstacles to enable the stem's unveiling from any obstacles between itself and the camera.

Figures <ref> - <ref> present the experimental results. In all figures, the green dotted line marks the start of the bimanual motion, that is at $7.67$ seconds. Figures <ref> - <ref> depict the successful reaching of the desired region around the stem as well as its unveiling with respect to surrounding occlusions. The camera's distance from the region's sphere $\vct{p}_r - \vct{p}_c$ decreases and becomes even less than $r$ (Fig. <ref>), that is the camera enters the ROI. Angle $\theta$ between $\vct{p}_r - \vct{p}_c$ and the camera’s view direction reaches a near-zero value (Fig. <ref>), thus the centering is achieved. Notice that just before the bimanual motion of the two arms starts (green dotted line), there is a discontinuity due to the recalculation of $\vct{p}_{sb}$, that is the center of ROI. Figure <ref> shows the visible stem points, that increase steadily until the start of the bimanual motion as the camera moves for unveiling the stem. The drop in the number of visible points during the bimanual motion is due to the fact that the grasping arm moves inside the camera's field of view, thus hiding the stem momentarily. Nevertheless, this is corrected by the camera arm's control  continuing the stem's unveiling reaching an even higher number of unveiled points at the end.

Reaching. Desired region radius $r = 0.35$ m - red dashed line.

Centering. Angle $\theta$ between $\vct{p}_r - \vct{p}_c$ and camera’s view direction $\vct{z}_c$.

Reaching with centering.

Unveiling. Visible stem points with respect to camera.

Position error.

Grasped crop/Arm's end-effector alignment with stretched stem.

Grasping arm's position error and orientation alignment.

Applied force. Desired force magnitude $f_d$ - red dashed line.

Figures <ref> - <ref> depict the grasping arm's related control objectives achievement. As seen in Fig. <ref>, the norm of the arm's end-effector position error projected in its corresponding action space $\vct{I}_{3 \times 3} - \vct{n}_c \vct{n}_c^\intercal$ as well as the angle between the arm’s end-effector/grasped crop’s direction $\vct{y}_g$ and the stretched stem’s $\vct{n}_c$ both become zero. Moreover, the stem is stretched successfully, since the desired force magnitude $f_d$ is reached (Fig. <ref>); notice that the spike at $7.67$ seconds (green dotted line) is caused by the transition from a static state to motion when the grasping arm's control is activated  (<ref>).

Figure <ref> shows the camera's viewpoint at the start (Fig. <ref>-(a)), at the start of the bimanual motion (Fig. <ref>-(b)) and at the end of the task (Fig. <ref>-(c)), enabling the stem cutting.

In-hand camera's viewpoint. (a) Process start. (b) Bimanual motion's start. (c) End of overall task.

§ CONCLUSIONS

In this work, a bimanual control methodology for reaching a pre-cut state regarding the target crop's stem is presented. Control methods are designed for end-effector velocities for the camera and the grasping arm in order to reach and unveil the stem  as well as manipulate the crop to generate space around the stem. Lab experimentation with a mock-up vine setup and a plastic grape cluster validates the methodology's success. Future work involves the application of the proposed method in a vineyard.

[1]
V. Marinoudi, C. G. Sørensen, S. Pearson, and D. Bochtis, “Robotics and
labour in agriculture. a context consideration,” Biosyst. Eng., vol.
184, pp. 111–121, Aug. 2019, doi:

[2]
M. Eberhardt and D. Vollrath, “The effect of agricultural technology on the
speed of development,” World Dev., vol. 109, pp. 483–496, Sept.
2018, doi: https://doi.org/10.1016/j.worlddev.2016.03.017.

[3]
F. A. Fathallah, “Musculoskeletal disorders in labor-intensive agriculture,”
Appl. Ergon., vol. 41, no. 6, pp. 738–743, 2010, doi:

[4]
S. Fountas, N. Mylonas, I. Malounas, E. Rodias, C. H. Santos, and E. Pekkeriet,
“Agricultural robotics for field operations,” Sensors, vol. 20,
no. 9, p. 2672, May 2020, doi: https://doi.org/10.3390/s20092672.

[5]
C. W. Bac, E. J. van Henten, J. Hemming, and Y. Edan, “Harvesting robots for
high-value crops: State-of-the-art review and challenges ahead,” Int.
J. Agric. Biol. Eng., vol. 31, no. 6, pp. 888–911, Nov. 2014, doi:

[6]
R. R. Shamshiri, C. Weltzien, I. A. Hameed, I. J. Yule, T. E. Grift, S. K.
Balasundram, L. Pitonakova, D. Ahmad, and G. Chowdhary, “Research and
development in agricultural robotics: A perspective of digital farming,”
Int. J. Agric. Biol. Eng., vol. 11, no. 4, pp. 1–11, 2018, doi:

[7]
Y. Onishi, T. Yoshida, H. Kurita, T. Fukao, H. Arihara, and A. Iwai, “An
automated fruit harvesting robot by using deep learning,” ROBOMECH
Journal, vol. 6, no. 13, pp. 2––9, 2019, doi:

[8]
A. Silwal, J. R. Davidson, M. Karkee, C. Mo, Q. Zhang, and K. Lewis, “Design,
integration, and field evaluation of a robotic apple harvester,”
Journal of Field Robotics, vol. 34, no. 6, pp. 1140––1159, Sept.
2017, doi: https://doi.org/10.1002/rob.21715.

[9]
S. Hayashi, K. Ganno, Y. Ishii, and I. Tanaka, “Robotic harvesting system for
eggplants,” Japan Agricultural Research Quarterly: JARQ, vol. 36,
no. 3, pp. 163––168, 2002, doi: https://doi.org/10.6090/jarq.36.163.

[10]
E. J. V. Henten, J. Hemming, B. van Tuijl, J. Kornet, J. Meuleman, J. Bontsema,
and E. van Os, “An autonomous robot for harvesting cucumbers in
greenhouses,” Autonomous Robots, vol. 13, pp. 241––258, 2002, doi:

[11]
Q. Feng, W. Zou, P. Fan, C. Zhang, and X. Wang, “Design and test of robotic
harvesting system for cherry tomato,” International Journal of
Agricultural and Biological Engineering, vol. 11, no. 1, pp. 96––100,
2018, doi: https://doi.org/10.25165/j.ijabe.20181101.2853.

[12]
S. Hayashi, K. Shigematsu, S. Yamamoto, K. Kobayashi, Y. Kohno, J. Kamata, and
M. Kurita, “Evaluation of a strawberry-harvesting robot in a field test,”
Biosystems Engineering, vol. 105, no. 2, pp. 160––171, Feb. 2010,
doi: https://doi.org/10.1016/j.biosystemseng.2009.09.011.

[13]
Y. Zhao, L. Gong, C. Liu, and Y. Huang, “Dual-arm robot design and testing for
harvesting tomato in greenhouse,” IFAC-PapersOnLine, vol. 49, no. 16,
pp. 161––165, 2016, doi: https://doi.org/10.1016/j.ifacol.2016.10.030.

[14]
X. Ling, Y. Zhao, L. Gong, C. Liu, and T. Wang, “Dual-arm cooperation and
implementing for robotic harvesting tomato using binocular vision,”
Rob. Auton. Syst., vol. 114, pp. 134––143, Apr. 2019, doi:

[15]
A. Sahbani, S. El-Khoury, and P. Bidaud, “An overview of 3d object grasp
synthesis algorithms,” Robotics and Autonomous Systems, vol. 60,
no. 3, pp. 326–336, 2012, doi: https://doi.org/10.1016/j.robot.2011.07.016.

[16]
O. Ringdahl, P. Kurtser, and Y. Edan, “Evaluation of approach strategies for
harvesting robots: Case study of sweet pepper harvesting: Category: (5),”
Journal of Intelligent & Robotic Systems, vol. 95, no. 1, pp.
149––164, July 2019, doi: https://doi.org/10.1007/s10846-018-0892-7.

[17]
B. Keinert, M. Innmann, M. Sänger, and M. Stamminger, “Spherical fibonacci
mapping,” ACM Transactions on Graphics, vol. 34, no. 6, pp. 1–7,
2015, doi: https://doi.org/10.1145/2816795.2818131.

[18]
C. Cheah and D. Wang, “Region reaching control of robots: Theory and
experiments,” in Proceedings of the 2005 IEEE International Conference
on Robotics and Automation, 2005, pp. 974–979.

[19]
T. Kastritsi, D. Papageorgiou, I. Sarantopoulos, S. Stavridis, Z. Doulgeri, and
G. A. Rovithakis, “Guaranteed active constraints enforcement on point
cloud-approximated regions for surgical applications,” in 2019
International Conference on Robotics and Automation (ICRA), 2019, pp.
8346–8352, doi: https://doi.org/10.1109/ICRA.2019.8793953.


# Deep Multiway Canonical Correlation Analysis for
Multi-subject EEG Normalization

###### Abstract

The normalization of brain recordings from multiple subjects responding to the
natural stimuli is one of the key challenges in auditory neuroscience. The
objective of this normalization is to transform the brain data in such a way
as to remove the inter-subject redundancies and to boost the component related
to the stimuli. In this paper, we propose a deep learning framework to improve
the correlation of electroencephalography (EEG) data recorded from multiple
subjects engaged in an audio listening task. The proposed model extends the
linear multi-way canonical correlation analysis (CCA) for audio-EEG analysis
using an auto-encoder network with a shared encoder layer. The model is
trained to optimize a combined loss involving correlation and reconstruction.
The experiments are performed on EEG data collected from subjects listening to
natural speech and music. In these experiments, we show that the proposed deep
multi-way CCA (DMCCA) based model significantly improves the correlations over
the linear multi-way CCA approach with absolute improvements of $0.08$ and
$0.29$ in terms of the Pearson correlation values for speech and music tasks
respectively.

Index Terms—  Canonical correlation analysis (CCA), multi-way CCA, Deep CCA,
Audio-EEG analysis.

##  1 Introduction

The analysis of brain signals as a response to an auditory stimulus is of
profound interest from both the understanding as well as the application
perspective. The main goal of this research is to find a relationship between
the stimulus and the response. The response signal in non-invasive
measurements is typically in the form of electroencephalography (EEG),
magnetoencephalography (MEG) or, functional magnetic resonance imaging (fMRI).
Among these modalities, the EEG represents the most easily accessible brain
recordings for analyzing speech comprehension and audio signal perception [1].

Though the EEG can provide high temporal resolution (typically sampled at
$128$ Hz or more), the data is captured at the scalp. As a result, the desired
response is significantly attenuated by the skull and makes the recordings
highly noisy (average SNR of $-20$dB [2]). For short stimuli, it is a common
practice to repeat the experiment for a given condition. The averaging
operation on the brain responses improves the SNR. This resultant average
response is termed as the event related potential (ERP) [3]. However, for
naturalistic auditory stimuli like speech and audio, repetition of the
experiments with the same conditions is often impractical. This has brought in
greater demand for single-trial EEG analysis methods.

The temporal response function (TRF) [4] is one of the earliest in this
direction of single-trial EEG analysis. This method assumes a linear
relationship between the stimuli and the EEG response. The system
identification is performed using least-squares and the performance of these
models is quantified by the Pearson correlation between the estimated signal
and the actual signal. In the EEG recordings, the presence of signals
unrelated to the presented stimuli, as well as the attenuation by the scalp,
results in shallow correlation values for these models in the range of
$0.1-0.2$ [4].

The canonical correlation analysis (CCA) tries to alleviate this problem to
some extent by finding common transformations of the stimuli and the
responses. It finds a pair of optimum linear transforms, for the stimulus and
the response, such that the correlation between the two new projections is
maximized [5]. In order to characterize the non-linear relationship between
the stimuli and EEG response, we had recently proposed a deep CCA [6] for
audio-EEG datasets [7].

Most of the single-trial analysis methods model the stimulus-response data
only for a single subject. One way to expand the dataset is to perform the
experiment on multiple subjects for the same stimuli and conditions. At the
same time, as many other subject-specific factors drive the EEG response, the
problem of finding common signals across the subjects’ responses is
challenging. The linear CCA can be performed only on two views of the data at
a time. A multi-variate version of the CCA, termed as multiway CCA or
generalized CCA [8, 9, 10], has been proposed for aggregating EEG response
from multiple-views (subjects). As all views (EEG responses) are based on the
same stimulus, the components that are common across the views [11] can elicit
the information regarding the brain processing of the stimuli. The multiway
CCA (MCCA) successfully extracts common signals across more than two data-
views [12]. However, all these models (TRF, CCA and MCCA) are based on linear
assumptions of the brain functions. In this paper, we hypothesize that a
linear system-response assumption may be too simplistic and inadequate to
uncover the transfer function of the brain responses.

In this paper, we propose a deep learning model for the multi-subject EEG
normalization task. Given the EEG response to the same stimuli from multiple
subjects, the proposed model projects the EEG responses to a space that
maximizes the inter-subject correlation. This model, termed as deep multiway-
CCA, extends the linear methods to allow the analysis of single-trial multi-
subject audio-EEG datasets. The model uses a reconstruction based
regularization to alleviate the issues of over-fitting in the deep MCCA model.
The deep MCCA model is compared with linear MCCA [12, 11] for speech [13] and
audio [14] perception tasks. In our experiments comparing the linear MCCA with
the proposed DMCCA, the DMCCA model improves correlation on the average by
$0.08$ for speech tasks and $0.18$ for audio tasks. These improvements in the
correlation values are found to be statistically significant ($p<0.01$) on a
pairwise t-test analysis as well.

##  2 Related Prior Work

The recent years have seen the use of deep learning for several brain mapping
tasks like computational memory prediction [15] and for reconstruction of
brain activity for visual stimuli [16]. A review of several efforts in
decoding brain activity using deep learning techniques is given in Zheng et
al. [17].

In auditory tasks, EEG recordings have shown to contain rhythm information in
music perception [18]. Recently, Das et al. [19] showed that deep learning
techniques can improve auditory attention decoding in the perception of noisy
speech. In multi-speaker cocktail party scenarios, Deckers et al. [20] showed
that neural networks can identify the attended speaker.

In the context of CCA, the extensions to connect multiple data matrices have
been proposed under names such as multiple CCA, multiway CCA, multiset CCA,
generalized CCA, or multiway orthonormalized partial least squares [21, 22,
23]. There have been attempts to generalize the deep CCA [24]. Liu et al. [25]
has proposed a deep model, referred to as DGCCA [24] in the literature, for
SSVEP frequency recognition. To the best of our knowledge, this paper presents
the first efforts in developing deep models for normalization of audio-EEG
data from multiple subjects.

##  3 Deep MCCA

###  3.1 Background - Linear Multi-way CCA

The multiway CCA (MCCA) is a linear method that expands the linear CCA to more
than two data-views. It obtains an optimum linear transform for each data-
view, such that the inter-set correlation (ISC) across all the projections is
maximum [11, 12].

Let $N$ data-views be denoted as $\mathbf{x}_{n}\in\mathcal{R}^{d_{n}}$ for
$n=1\text{ to }N$ and $D_{N}=\sum_{n}d_{n}$. The MCCA tries to extract the
signals common across the data-views. For $1$D projection, the MCCA projects
the $N$ data-views onto a $1$D subspace, such that all the views are highly
correlated to each other. Let, the vector
$\mathbf{v}_{n}\in\mathcal{R}^{d_{n}\times 1}$ denote projection vector that
linearly transforms data-view $\mathbf{x}_{n}$ onto the $1$D subspace. The
linear transforms are learnt such that the inter-set correlation (ISC) among
the projections is maximum. The ISC can be formulated as:

$\text{ISC}=\frac{1}{N-1}\frac{r_{B}}{r_{W}}$ (1)

where $r_{B}$ and $r_{W}$ are the between-set covariance and within-set
covariance respectively given by,

$\displaystyle r_{B}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1,j\neq
i}^{N}\mathbf{v}_{i}^{\top}\mathbf{R}^{ij}\mathbf{v}_{j}$ (2) $\displaystyle
r_{W}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\mathbf{v}_{i}^{\top}\mathbf{R}^{ii}\mathbf{v}_{i}.$
(3)

Here, $\mathbf{R}^{ij}\in\mathcal{R}^{d_{i}\times d_{j}}$ are the cross-
correlation matrix of the data-views $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$.

The optimal linear transform vectors $\\{\mathbf{v}_{n}\\}_{n=1}^{N}$ are
obtained by finding the eigenvector with the maximum eigenvalue of the
following eigen problem [11]:

$\mathbf{R}\mathbf{v}=\lambda\mathbf{D}\mathbf{v}$ (4)

where $\lambda=(N-1)\rho+1$, $\mathbf{v}\in\mathcal{R}^{D_{N}\times 1}$ is the
concatenated vector of the linear transform vectors
$\\{\mathbf{v}_{n}\\}_{n=1}^{N}$, matrix
$\mathbf{R}\in\mathcal{R}^{D_{N}\times D_{N}}$ is a block matrix with its
block elements as $[\mathbf{R}]_{ij}=\mathbf{R}^{ij}$ and
$\mathbf{D}\in\mathcal{R}^{D_{N}\times D_{N}}$ is a block-diagonal matrix with
its elements as $[\mathbf{D}]_{ii}=\mathbf{R}^{ii}$. If the $N$ data-views are
to be projected onto a $d$ dimensional ($d>1$) subspace, then the linear
transform matrix for each data-view $\mathbf{x}_{n}$ is
$\mathbf{V}_{n}\in\mathcal{R}^{d_{n}\times d}$, obtained by choosing the top
‘$d$’ eigenvectors of Eq. 4 .

###  3.2 Deep Multiway CCA

####  3.2.1 Natural Extension - DGCCA

Using the $N$ random multi-variates, the goal of DGCCA is to determine $N$
neural transformations (non-linear) such that the correlation cost is
maximized. Let $f_{n}(\cdot)$, for $n=1\text{ to }N$, denote the neural
network that takes $\mathbf{x}_{n}$ as input and outputs a $d$ dimension
representation. Let each neural network’s trainable parameters be denoted as
$\boldsymbol{\theta}_{n}$. The parameters are learned such that the sum of all
inter-variate correlations is maximized. This correlation cost is,

$\displaystyle\rho=\sum_{i=1}^{N}\sum_{j=1,j\neq
i}^{N}\operatorname{corr}\left(f_{i}\left(\mathbf{x}_{i};\boldsymbol{\theta}_{i}\right),f_{j}\left(\mathbf{x}_{j};\boldsymbol{\theta}_{j}\right)\right)$
(5)

We refer to this model as the deep generalized canonical correlation analysis
(DGCCA) model in Figure 1.

Fig. 1: The DGCCA and the proposed DMCCA.

####  3.2.2 Proposed Model

The $N$ neural networks in the DGCCA model transform the $d_{n}$ dimensional
data-views to $d$ ($d\leq d_{n}$) dimensional common subspace. This
transformation can be viewed as a shared encoder of the views. In the proposed
model, we also construct a decoder which attempts to reconstruct the data-
views from the shared encoder output. This architecture is also illustrated in
Figure 1.

The encoder model consists of $N$ neural transformations, each operating on
one of the $N$ data-views. The $N$ data-views, propagated through the encoding
parts, generate $N$ vectors of $d$ dimensions
($f_{i}\left(\mathbf{x}_{i};\boldsymbol{\theta}_{i}\right)$) which are
concatenated to form a $Nd$ dimensional vector (shared encoder output). This
$Nd$ dimensional encoder output ($\mathbf{y}$) is used as the common
representation for all the decoder layers. The decoder also has $N$ parallel
networks, each of which attempts to reconstruct the corresponding data-view
$\mathbf{x}_{n}$. This proposed model is referred to as the deep multiway
canonical correlation analysis (DMCCA) model.

The DMCCA model is trained to maximize the inter-set correlation among the
encoders’ outputs, and minimize the MSE (Mean Square Error) loss at each of
decoder outputs. Let the trainable parameters of the model be split into $N$
parts, where each
$\\{\boldsymbol{\theta}_{n}\,,\boldsymbol{\hat{\theta}}_{n}\\}_{n=1}^{N}$
corresponds to the parameters of the encoder and decoder of the
$n^{\text{th}}$ data-view. Then, the optimization of the DMCCA model is
performed by maximizing

$\displaystyle\rho^{\prime}=\rho-\lambda\sum_{n=1}^{N}\operatorname{MSE}\left(\mathbf{x}_{n},\hat{f}_{n}(\mathbf{y};\boldsymbol{\hat{\theta}}_{n})\right)$
(6)

where $\rho$ is defined by Eq. (5) and $\operatorname{MSE}(\cdot)$ denotes the
mean squared error between the reconstructed data-view and the input data-view
(auto-encoding loss). The parameter $\lambda$ is a hyper-parameter that
controls the relative importance of the two-cost functions. It is found that
the MSE regularization loss is complementary to the inter-variate correlation
loss. Further, the DGCCA is a special case of DMCCA for $\lambda=0$.

Following the model training, the encoder part of each data-view is used as
the fixed non-linear transform that normalizes each data-view to remove
subject-specific redundancies in the data and boost the stimuli related
components.

##  4 Audio-EEG Setup

###  4.1 Datasets

We perform experiments on two naturalistic auditory stimuli and their
respective EEG responses: a speech dataset and an audio dataset. The speech
dataset is taken from a study on cortical entrainment to natural speech [13].
The audio dataset is the NMED-H which consists of EEG responses from $48$
subjects listening to natural music [14].

The speech dataset is the same dataset used for the baseline linear MCCA model
[12]. This dataset is also used in the application of CCA [26] and deep CCA
[7] models. The EEG data are recorded using $128$ electrodes. We perform the
same pre-processing steps done by Cheveigné et al. [12]. The EEG data are
down-sampled to $64$ Hz and de-trended and de-noised using noise tools
software [27]. Then, EEG data are passed through a band-pass filter with a
pass-band in the range of $0.1-12$ Hz. The stimuli data are obtained from
audio sampled at $44100$ Hz. The audio envelope is squared and smoothed by
convolving the envelope with a square window [27]. Finally, the stimuli data
(audio envelopes) are downsampled to $64$ Hz and are further compressed with a
power of $1/3$. More details about the pre-processing can be found in
Cheveigné et al. [26].

The audio dataset is the NMED-H dataset [14]. It is designed to study the
brain signals for a natural music listening task. It contains EEG responses to
intact and scrambled versions of $4$ full-length Hindi pop songs. Each
stimulus is approximately $4.5$ minutes long, sampled at $44.1k$Hz. The EEG
data are recorded using $125$ electrodes. They are sampled at $125$ Hz with
average reference. More details of the data acquisition and pre-processing are
given in Kaneshiro [28]. The features from the audio stimuli are extracted in
a similar way as described by Gang et al. [14]. The acoustic features are
extracted using the music information retrieval (MIR) toolbox [29]. We have
extracted $20$ dimensional features in $25$ms analysis windows with a $50$%
overlap between frames [30]. The Principal component analysis (PCA) is
performed on these $20$D data to extract a $1$D representation (PC$1$) for the
data. Along with the PC$1$, two individual features, root mean square (RMS)
and spectral flux, that reflect amplitude envelope and timbre respectively are
chosen to obtain a $3$D representation for the stimuli. The final sampling
rate of the acoustic features is $80$ Hz. Similarly, the EEG response is
resampled to $80$ Hz. All experiments on the audio dataset are performed on
each of these $3$D pre-processed features or the $1$D envelope features.

###  4.2 Performance metric

The outputs of linear/deep MCCA models are provided to a linear CCA based
model. It is applied to further transform the MCCA/DMCCA model outputs and the
stimulus data, separately for each subject. We use the CCA3 configuration
described in Cheveigné et al. [26].

The performance comparison between the linear and deep versions of the MCCA
models is done using the final correlation scores of the first canonical
covariate of the CCA3 method. Another metric to measure the models’
performance is the classification accuracy. For that, we have performed a
Cohen’s d-prime match-vs-mismatch classification analysis as discussed by
Cheveigné et al. [26].

###  4.3 DMCCA Model

For EEG response from $N$ subjects and a $d_{s}$ time-lagged $1$D (common)
stimulus, the encoder section consists of two hidden layers with $60$ nodes
each and an output layer of $10$ nodes. Thus, the final de-noised
representations for the EEG responses is $10$D. The $N+1$ ($N$ EEG responses
plus the stimuli features) vectors of $10$D are concatenated to form a
$10(N+1)$D vector. This vector is provided to each decoder. The decoder
section also consists of two hidden layers of $60$ nodes and $110$ nodes
respectively.

The final outputs from the linear MCCA are of $128$D and $125$D for speech and
audio data respectively (same as the input). The outputs from the deep MCCA
are of $10$D. The linear MCCA first encodes the data onto $d_{s}$ dimensions
and then projects them onto the $128$D space [12] (for speech dataset) or
$125$D  (  for the audio dataset  ) . The choice of time-lag ($d_{s}$) for
each version of MCCA is done by varying it from $10$ to $110$. The time-lag
which provided the best performance is chosen.

To make the final $1$D representations of both the MCCA versions comparable,
the $10$D denoised EEG response of each subject is provided to the filterbank
from the CCA$3$ method [26] and a PCA to project them onto a $139$D subspace.
The $10$D stimulus features are projected onto a $1$D subspace using PCA and
provided to the filterbank, resulting in a $21$D stimulus features. These
features are provided to a linear CCA model to obtain $1$D representations
[26].

###  4.4 Experimental Setup

Fig. 2: The correlation values for linear and deep MCCA, on $6$ subjects from
the speech dataset. Pairwise t-test statistical significance is shown (ns
implies no significance ($p>0.01$), ** implies $p\leq 0.01$), *** implies
$p\leq 0.001$), **** implies $p\leq 1e-4$)

From the speech dataset, stimulus-response of $6$ subjects is used in the
experiments. Each subject has $20$ recorded sessions with approximately $160$
seconds in each session. We perform $20$ fold-validation experiments on these
$20$ sessions. In each validation experiment, $18$ of the $20$ sessions are
used for training the models, while one of the remaining sessions is used for
validation and the other for testing. This results in about $200k$ samples for
training the model.

The audio dataset has EEG response from $48$ subjects listening to $4$
versions of $4$ full length Hindi pop songs, each having about $4.5$ minutes
duration. We aggregated the stimulus-response data based on the common
stimuli. Each subject has $2$ EEG recording trials for $2$ different stimuli.
This results in about $38k$ samples for training and $2k$ samples for
validation and testing respectively.

The $N$ subjects EEG responses and a time-lagged stimulus (time lag of $40$
for linear and $60$ for deep MCCA) features are provided to the linear and
deep versions of MCCA. The denoised outputs of the MCCA (linear/deep) are
provided to the CCA3 method [26]. The final representations of $1$D are
obtained by the CCA$3$ method and their correlations are used for comparing
the MCCA methods.

A leaky ReLU activation function with a negative slope coefficient of $0.1$ is
used at the encoder layer in DMCCA model. A linear activation function is used
at the output layer of the decoder sections. The dropout regularization [31]
is incorporated in the deep models, and the dropout with the best perfomance
is chosen.

The MCCA (linear/deep) denoising step is performed on the combination of EEG
responses whereas the final step of CCA3 method is done for each subject
separately.

##  5 Results

###  5.1 Performance analysis

For speech dataset, the linear and deep MCCA methods are performed on a
collection of $6$ subjects. The results of the cross-validation experiments
for all the $6$ subjects are presented in Fig. 2. As seen in the Fig. 2, the
audio-EEG correlation is improved for all the subjects using DMCCA. On the
average, the DMCCA improves the correlation by around $0.08$. Additional
experiments comparing the performance of the DGCCA model [24], proposed by Liu
et al. [25] and our proposed DMCCA model has been done on the speech dataset.
The average correlation of the the six subjects has increased from $0.158$
(for the DGCCA model) to $0.234$ (our proposed model). This shows that this
MSE regularization term helps the model improve the correlations.

For the audio dataset, the $48$ subjects were divided into sets based on the
common stimulus, with $12$ subjects in each set. The MCCA methods were tried
on these sets separately. One experiment with the preprocessed envelope
representing the stimuli features and $3$ experiments with the $3$D features
(from 4.1) were performed for both versions of MCCA. The results are presented
in the Table. 1. The DMCCA consistently improves over linear MCCA, on an
average of $0.29$ in the correlation values.

Table 1: Average correlation values of $48$ subjects from NMED-H Dataset | LMCCA | DMCCA
---|---|---
Envelope | 0.0694 | 0.3268
PC1 | 0.0237 | 0.3552
RMS | 0.0305 | 0.3169
Spectral Flux | 0.0432 | 0.3261
Table 2: The d-prime metric for both the datasets (in percentages) | Speech dataset | NMED-H
---|---|---
segment size | LMCCA | DMCCA | LMCCA | DMCCA
1s | 2.84 | 2.91 | 1.36 | 2.12
2s | 3.02 | 4.27 | 3.45 | 4.80
4s | 5.39 | 7.52 | 7.18 | 11.84
8s | 8.45 | 11.29 | 14.94 | 34.18
16s | 15.87 | 18.57 | 32.38 | 38.39
32s | 19.12 | 32.90 | - | -

###  5.2 Statistical Analysis

To test the statistical significance, one tail pairwise t-test is performed on
the final correlation values for the linear and deep versions of MCCA. For the
speech dataset, it is seen that the improvements are statistically significant
($p<1e-2$) for $5$ out of $6$ subjects. It is found to be statistically
significant ($p<1e-2)$ for all the $4$ stimulus features for the audio
dataset. Similar to Cheveigné et al. [26], Cohen’s d-prime value for match-vs-
mismatch classification analysis is also performed. Segments of EEG and
stimuli of duration $T$ are selected and the correlation among all the
combinations of EEG and stimuli are measured. The correlation of an EEG
response and its related stimulus falls under the ’match’ class, whereas for
unrelated stimuli fall under ’mismatch’ class. For $T$ varying from $1$ to
$32$ seconds, the d-prime values for the MCCA methods, for both the datasets,
are tabulated in Table. 2. As seen in Table. 2, for all the segment durations,
the deep MCCA improves over the linear version for both the datasets.

##  6 Summary

In this paper, we have proposed a deep multi-way CCA framework for multi-
subject EEG analysis of auditory stimuli. The DMCCA model uses the EEG data
from multiple subjects in an autoencoder setup with a shared layer. The model
uses a combination of reconstruction loss and correlation loss. The
experiments on the speech dataset shows that the model significantly improves
the correlations over the linear multi-way CCA. A classification based measure
of stimuli-response alignment also shows the benefits of the proposed DMCCA
model. Further, we show that the improvements seen in speech data are
consistent for EEG-music datasets as well.

## References

  * [1] Malte Wöstmann, Lorenz Fiedler, and Jonas Obleser,  “Tracking the signal, cracking the code: Speech and speech comprehension in non-invasive human electrophysiology,”  Language, Cognition and Neuroscience, vol. 32, no. 7, pp. 855–869, 2017.
  * [2] Saeid Sanei and Jonathon A Chambers,  “EEG signal processing,”  Wiley Online Library, 2007.
  * [3] Michael GH Coles and Michael D Rugg,  Event-related brain potentials: An introduction.,  Oxford University Press, 1995.
  * [4] Edmund C Lalor, Alan J Power, Richard B Reilly, and John J Foxe,  “Resolving precise temporal processing properties of the auditory system using continuous stimuli,”  Journal of neurophysiology, vol. 102, no. 1, pp. 349–359, 2009\.
  * [5] Bruce Thompson,  Canonical correlation analysis: Uses and interpretation,  Number 47. Sage, 1984.
  * [6] Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu,  “Deep canonical correlation analysis,”  in International Conference on Machine Learning, 2013, pp. 1247–1255.
  * [7] Jaswanth Reddy Katthi, Sriram Ganapathy, Sandeep Kothinti, and Malcolm Slaney,  “Deep canonical correlation analysis for decoding the auditory brain,”  in 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE, 2020, pp. 3505–3508.
  * [8] Nicolle M Correa, Tom Eichele, Tülay Adalı, Yi-Ou Li, and Vince D Calhoun,  “Multi-set canonical correlation analysis for the fusion of concurrent single trial erp and functional mri,”  Neuroimage, vol. 50, no. 4, pp. 1438–1445, 2010.
  * [9] Xiao Fu, Kejun Huang, Mingyi Hong, Nicholas D Sidiropoulos, and Anthony Man-Cho So,  “Scalable and flexible multiview max-var canonical correlation analysis,”  IEEE Transactions on Signal Processing, vol. 65, no. 16, pp. 4150–4165, 2017.
  * [10] Qiong Zhang, Jelmer P Borst, Robert E Kass, and John R Anderson,  “Inter-subject alignment of meg datasets in a common representational space,”  Human brain mapping, vol. 38, no. 9, pp. 4287–4301, 2017.
  * [11] Lucas C Parra,  “Multi-set canonical correlation analysis simply explained,”  arXiv preprint arXiv:1802.03759, 2018.
  * [12] Alain de Cheveigné, Giovanni M Di Liberto, Dorothée Arzounian, Daniel DE Wong, Jens Hjortkjær, Søren Fuglsang, and Lucas C Parra,  “Multiway canonical correlation analysis of brain data,”  NeuroImage, vol. 186, pp. 728–740, 2019.
  * [13] Giovanni M Di Liberto, James A O’Sullivan, and Edmund C Lalor,  “Low-frequency cortical entrainment to speech reflects phoneme-level processing,”  Current Biology, vol. 25, no. 19, pp. 2457–2465, 2015.
  * [14] Nick Gang, Blair Kaneshiro, Jonathan Berger, and Jacek P Dmochowski,  “Decoding neurally relevant musical features using canonical correlation analysis.,”  in ISMIR, 2017, pp. 131–138.
  * [15] Xuyun Sun, Cunle Qian, Zhongqin Chen, Zhaohui Wu, Benyan Luo, and Gang Pan,  “Remembered or forgotten?—an eeg-based computational prediction approach,”  PloS one, vol. 11, no. 12, pp. e0167497, 2016.
  * [16] Ye Yuan, Guangxu Xun, Qiuling Suo, Kebin Jia, and Aidong Zhang,  “Wave2vec: Deep representation learning for clinical temporal data,”  Neurocomputing, vol. 324, pp. 31–42, 2019.
  * [17] Xiao Zheng, Wanzhong Chen, Mingyang Li, Tao Zhang, Yang You, and Yun Jiang,  “Decoding human brain activity with deep learning,”  Biomedical Signal Processing and Control, vol. 56, pp. 101730, 2020\.
  * [18] Sebastian Stober, Daniel J Cameron, and Jessica A Grahn,  “Using convolutional neural networks to recognize rhythm stimuli from electroencephalography recordings,”  in Advances in neural information processing systems, 2014, pp. 1449–1457.
  * [19] Neetha Das, Jeroen Zegers, Tom Francart, Alexander Bertrand, et al.,  “Linear versus deep learning methods for noisy speech separation for eeg-informed attention decoding,”  Journal of Neural Engineering, vol. 17, no. 4, pp. 046039, 2020\.
  * [20] Lucas Deckers, Neetha Das, Amir Hossein Ansari, Alexander Bertrand, and Tom Francart,  “Eeg-based detection of the attended speaker and the locus of auditory attention with convolutional neural networks,”  bioRxiv, p. 475673, 2018.
  * [21] Samuel M Gross and Robert Tibshirani,  “Collaborative regression,”  Biostatistics, vol. 16, no. 2, pp. 326–338, 2015.
  * [22] Irene Sturm,  Analyzing the perception of natural music with EEG and ECoG,  Ph.D. thesis, Berlin, Technische Universität Berlin, 2016.
  * [23] Yu Zhang, Guoxu Zhou, Qibin Zhao, Akinari Onishi, Jing Jin, Xingyu Wang, and Andrzej Cichocki,  “Multiway canonical correlation analysis for frequency components recognition in ssvep-based bcis,”  in International conference on neural information processing. Springer, 2011, pp. 287–295.
  * [24] Adrian Benton, Huda Khayrallah, Biman Gujral, Dee Ann Reisinger, Sheng Zhang, and Raman Arora,  “Deep generalized canonical correlation analysis,”  arXiv preprint arXiv:1702.02519, 2017.
  * [25] Qianqian Liu, Yong Jiao, Yangyang Miao, Cili Zuo, Xingyu Wang, Andrzej Cichocki, and Jing Jin,  “Efficient representations of eeg signals for ssvep frequency recognition based on deep multiset cca,”  Neurocomputing, vol. 378, pp. 36–44, 2020.
  * [26] Alain de Cheveigné, Daniel DE Wong, Giovanni M Di Liberto, Jens Hjortkjaer, Malcolm Slaney, and Edmund Lalor,  “Decoding the auditory brain with canonical component analysis,”  NeuroImage, vol. 172, pp. 206–216, 2018.
  * [27] Alain de Cheveigné and Dorothée Arzounian,  “Robust detrending, rereferencing, outlier detection, and inpainting for multichannel data,”  NeuroImage, vol. 172, pp. 903–912, 2018.
  * [28] Blair Bohannan Kaneshiro,  Toward an Objective Neurophysiological Measure of Musical Engagement,  Ph.D. thesis, Stanford University, 2016.
  * [29] Olivier Lartillot and Petri Toiviainen,  “A matlab toolbox for musical feature extraction from audio,”  in International conference on digital audio effects. Bordeaux, 2007, pp. 237–244.
  * [30] George Tzanetakis and Perry Cook,  “Musical genre classification of audio signals,”  IEEE Transactions on speech and audio processing, vol. 10, no. 5, pp. 293–302, 2002.
  * [31] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov,  “Dropout: a simple way to prevent neural networks from overfitting,”  The journal of machine learning research, vol. 15, no. 1, pp. 1929–1958, 2014.