# Multi-Modal Model Predictive Control through Batch Non-Holonomic Trajectory
Optimization: Application to Highway Driving

Vivek K. Adajania, Aditya Sharma, Anish Gupta, Houman Masnavi, K Madhava
Krishna and Arun K.Singh Manuscript received: September, 9, 2021; Revised
December, 6, 2022; Accepted January, 13, 2022.This paper was recommended for
publication by Stephen J. Guy upon evaluation of the Associate Editor and
Reviewers’ comments. This work was supported in part by the European Social
Fund via IT Academy program in Estonia, smart specialization project with BOLT
and grants COVSG24 and PSG605 from Estonian Research Council.The first three
authors are with the Robotics Research Center, Kohli Center for Intelligent
Systems, IIIT Hyderabad, India and Houman Masnavi and Arun Singh are with the
Institute of Technology, University of Tartu.Digital Object Identifier (DOI):
see top of this page.

###### Abstract

Standard Model Predictive Control (MPC) or trajectory optimization approaches
perform only a local search to solve a complex non-convex optimization
problem. As a result, they cannot capture the multi-modal characteristic of
human driving. A global optimizer can be a potential solution but is
computationally intractable in a real-time setting. In this paper, we present
a real-time MPC capable of searching over different driving modalities. Our
basic idea is simple: we run several goal-directed parallel trajectory
optimizations and score the resulting trajectories based on user-defined meta
cost functions. This allows us to perform a search over several locally
optimal motion plans. Although conceptually straightforward, realizing this
idea in real-time with existing optimizers is highly challenging from
technical and computational standpoints. With this motivation, we present a
novel batch non-holonomic trajectory optimization whose underlying matrix
algebra is easily parallelizable across problem instances and reduces to
computing large batch matrix-vector products. This structure, in turn, is
achieved by deriving a linearization-free multi-convex reformulation of the
non-holonomic kinematics and collision avoidance constraints. We extensively
validate our approach using both synthetic and real data sets (NGSIM) of
traffic scenarios. We highlight how our algorithm automatically takes lane-
change and overtaking decisions based on the defined meta cost function. Our
batch optimizer achieves trajectories with lower meta cost, up to 6x faster
than competing baselines.

##  I Introduction

Human driving is a complex mixture of discrete level decisions (merge,
overtake, etc.) and lower-level motion commands [comb_ad_planning]. If we
adopt an optimization perspective, the multiple discrete decisions can be seen
as local minima associated with the underlying non-convex trajectory
optimization problem [comb_ad_planning] [ad_logical_const]. Local optimizers
based on Sequential Quadratic Programming (SQP) or Gradient Descent (GD) are
not equipped to search over all the local minima. On the other hand, global
optimization techniques like mixed-integer programming [ad_logical_const]
offer a potential solution but are not particularly useful in a real-time
setting, especially in dense traffic scenarios.

Figure 1: The different colored samples represent the several locally optimal
trajectories for driving as close as possible to the cruise speed (meta cost).
The ego-vehicle is shown in red while the neighboring obstacles are shown in
blue. The number within the ellipses represent the velocities of the
respective entity. Unlike [fernet_planner], [behavior_ad], we explicitly
consider collision avoidance and kinematic constraints while generating
candidate trajectories. The trajectory shown in bold red achieves the best
performance (lowest meta-cost). We recommend seeing the accompanying video
(https://tinyurl.com/3wew7vu7) before reading the paper.

Main Idea: Let us define goals as tuples of position, velocity and
acceleration to be achieved by the ego-vehicle at the end of the planning
horizon. Then, our approach in this paper is built on a simple insight that
many different goal-directed trajectories can accomplish a given high-level
driving task. For example, if the task is to drive close to maximum velocity,
the autonomous car can choose its next goal to be in any lane. Some particular
goal choices may require overtaking a slow-moving vehicle directly in front,
while some may require safely merging with oncoming cars in a different lane.
More concretely, each goal-directed trajectories may converge to a different
local minima resulting in a multi-modal driving behavior (see pp-5
[comb_ad_planning]). Thus our proposed work is based on the idea of running
several parallel goal-directed trajectory optimization problems and ranking
the resulting locally optimal trajectories based on a user-defined meta cost
function. Although conceptually simple, we are not aware of any such approach
in existing works. There are trajectory sampling approaches such as
[fernet_planner], [behavior_ad] but they ignore collision avoidance and
kinematic constraints while computing possible maneuvers. One possible reason
the parallel/batch trajectory optimization approach has not been tried is that
solving non-convex problem associated with autonomous driving is challenging.
Running several instances of the problem in parallel only further increases
the computational burden. The conceptually straightforward approach of running
different optimizations in parallel CPU threads is not scalable for a large
batch size in a dense highway driving scenario (see Fig.2 and discussions
around it). We present a possible solution whose main novelties can be
summarized as follows.

Algorithmic: We present the very first batch non-holonomic trajectory
optimizer for real-time generation of several goal-directed locally optimal
collision-free trajectories in parallel. The core algorithmic challenge lies
in achieving linear scalability with respect to the batch size. As shown in
Fig.2(b), we cannot achieve such scalability by simply running each
optimization problem in a separate CPU thread. Instead, we need to parallelize
the batch optimizer’s per-iteration computation efficiently. We address the
scalability issue by developing a batch optimizer wherein iterating over the
different problem instances in parallel boils down to multiplying a single
constant matrix with a set of vectors. We show that the heavily vectorized
structure of our optimizer stems naturally from two key algorithmic
developments. First, we adopt a linearization-free multi-convex reformulation
of the kinematic and collision avoidance constraints. Second, we apply the
Alternating Minimization (AM) technique to solve the reformulated problem. We
rank the output of batch optimization with some simple meta cost functions
that model the higher-level driving objectives (e.g., driving with high-
speed). We show that the ranking mechanism coupled with an intelligent goal-
sampling approach automatically leads to discovering lane-change, vehicle
following, overtaking maneuvers based on the traffic scenario. Refer Section
II-B for a summary of algorithmic contribution over the author’s prior work.

Applied: We provide an open source implementation [github_repo_driving] for
review and to promote further research in this direction.

State-of-the-art Performance: We compare our batch optimizer based MPC with
three strong baselines and show that we achieve better solutions (in terms of
meta-cost value) while being up to 6x faster. Our first baseline is a standard
MPC that computes just a single locally optimal trajectory. Our second
baseline is batch multi-threaded implementation of optimal control solver
ACADO [acado]. Our final baseline is based on trajectory sampling in the
Frenet frame [fernet_planner].

Table I: Important Symbols  $x_{l}(t),y_{l}(t),\psi_{l}(t)$ |  Position and heading of the ego-vehicle at time $t$.
---|---
$\xi_{xj}(t),\xi_{yj}(t)$ |  Position of the $j^{th}$ obstacle at time $t$
$\alpha_{j,l}(t),d_{j,l}(t)$ |  Variables associated with our collision avoidance model. Refer to text for details.

##  II Preliminaries and Related Work

###  II-A Symbols and Notations

Small-case normal and bold font letters will be used to denote scalars and
vectors respectively. Bold-font upper-case letters will represent matrices.
The superscript $T$ will denote transpose of a matrix and vector. Some of the
main symbols are summarized in Table I while some are also defined in their
first place of use.

###  II-B Batch Non-Holonomic Trajectory Optimization

We are interested in solving $l$ non-holonomic trajectory optimizations in
parallel each of which can be formulated in the following manner. The
resulting trajectories from the parallel problems will be ranked based on a
meta cost function discussed later.

$\displaystyle\min\sum_{t}{\ddot{x}_{l}(t)^{2}+\ddot{y}_{l}(t)^{2}+\ddot{\psi}_{l}(t)^{2}}$
(1a)
$\displaystyle\dot{x}_{l}(t)=v_{l}(t)\cos\psi_{l}(t),\dot{y}_{l}(t)=v_{l}(t)\sin\psi_{l}(t),\forall
t$ (1b) $\displaystyle(x_{l}(t),y_{l}(t),\psi_{l}(t))\in\mathcal{C}_{b,l}$
(1c) $\displaystyle v_{min}\leq v_{l}(t)\leq
v_{max},\sqrt{\ddot{x}_{l}(t)^{2}+\ddot{y}_{l}(t)^{2}}\leq a_{max}$ (1d)
$\displaystyle-\frac{(x_{l}(t)-\xi_{xj}(t))^{2}}{a^{2}}-\frac{(y_{l}(t)-\xi_{yj}(t))^{2}}{b^{2}}+1\leq
0,$ (1e)

The subscript $l$ denotes that the specific variable belongs to the $l^{th}$
instance of the problem in the batch. The variables of the trajectory
optimization are $(x_{l}(t),y_{l}(t),\psi_{l}(t),v_{l}(t))$. The cost function
minimizes the squared acceleration value for the linear and angular motions.
The equality constraints (1b) stems from the non-holonomic kinematics of the
car. Constraints (1c) ensures the boundary conditions on the position, heading
angle and their derivatives. The inequalities (1d) represent the bounds on the
forward velocities and total acceleration. The set of constraints (1e)
enforces the collision avoidance between the ego and the neighboring vehicles
with the assumption that both are represented by road center-line aligned
ellipses. For the ease of exposition and without loss of generality, we assume
that every obstacle ellipse has the same major ($a$) and minor axis ($b$)
dimension. The ($a,b$) includes the inflation to account for the size of the
ego-vehicle. The ellipse of the ego-vehicle and obstacle will not overlap as
long the minimum separation distance is greater than $\sqrt{a^{2}+b^{2}}$. It
is worth pointing out that $a,b$ will be larger than the length and width of
the cars. In other words, there will be some over-approximation of the ego-
vehicle and obstacle footprints. To keep this value limited, we enforce some
restriction on the heading of the ego-vehicle. We discuss this more in the
beginning of Section III.

The typical control inputs for the ego-vehicle are the acceleration and the
steering inputs. The former can be obtained by the derivative or finite
difference of $v(t)$. The steering angle is given by
$\arctan(\frac{\dot{\psi}(t)*h}{v(t)})$, where $h$ is the inter-axle distance
[alonso_mora_parallel].

Note: All $l$ instances of (1a)-(1e) have the same velocity and acceleration
bounds and neighboring obstacles.

Existing Works: Trajectory optimizations of the form (1a)-(1e) are typically
solved through approaches like SQP [alonso_mora_parallel], GD [gd_ad] etc.
Authors in [ad_logical_const] proposed a global optimization approach but
considered the ego vehicle as a holonomic triple-integrator system. They
derived some approximations for the non-holonomic constraints that hold more
naturally at high forward velocities. In contrast, [knoll_ad] adopts a more
rigorous approach but the resulting algorithm was tested on environments
sparsely filled with obstacles. Authors in [2stage_opt] used an approach
similar to [ad_logical_const] to compute the best driving modality and then
refined the solution further through local optimization (e.g SQP) to handle
kinematic constraints.

Improvements over our prior efforts: Our approach extends recent work
[aks_iros21] to batch setting and further applies it to highway driving.
Specifically, we propose the core batch solution update rule (Eqn.(22), (24) )
by leveraging the implicit structures in the matrix algebra of [aks_iros21].
Furthermore, unlike our current work, [aks_iros21] did not consider
acceleration bounds. Typically, constraints on acceleration are modeled as
affine inequalities. However, such representation is not suitable for our
formulation that relies on reducing the trajectory optimization to a sequence
of unconstrained QPs to obtain an efficient batch update rule. We thus
reformulate the acceleration bounds in the same form as the collision
avoidance constraints of [aks_iros21]. Our current work is also an improvement
over [aks_icra20] that handled collision avoidance constraints through a novel
linearization approach. But as discussed in Section III-A, any linearization
of the underlying costs and constraints substantially complicates the
parallelization of the batch optimization.

##  III Main Results

This section derives our main algorithmic results. We begin by summarizing our
main assumptions.

1.Road Attached Frame: We assume that our problem set-up (i.e optimization
(1a)-(1e)) is defined in the reference frame of the center-line of the road
[2stage_opt]. This allows us to essentially treat curved roads as ones with a
straight-line geometry. The non-holonomic constraints are defined in the road
attached frame and holds true when the trajectories are reverted back to
global frame [2stage_opt].

2\. Restricted Heading Change: We post-process and discard locally optimal
trajectories resulting from (1a)-(1e) that incur a large heading change
($\approx\pm 13deg.$) with respect to the road center-line. The said heading
restriction is realistic in the case of typical highway driving scenarios
[nikolce_ad] and is made to accommodate the collision avoidance model defined
in (1e).

3\. Trajectory Prediction: The batch optimization generates trajectories that
are feasible with respect to the constant-velocity prediction of the
trajectories of the dynamic obstacles (neighboring cars). That is,
$\xi_{x,j}(t),\xi_{y,j}(t)$ are obtained by linearly interpolating the
positions with the current velocity. We adopt such a minimalist representation
to test the full potential of our batch optimization, specifically how its
fast re-planning ensures safety in the absence of any complex trajectory
forecasting algorithms.

###  III-A Intuition from Parallel Least Squares

Consider the following $l$ linear least squares problem

$\displaystyle\min_{\textbf{s}_{l}}\|\textbf{F}\textbf{s}_{l}-\textbf{g}_{l}\|_{2},\forall
l=1,2....l$ (2)

All problem instances share the same matrix F but have different vector
$\textbf{g}_{l}$. A simple way to solve the problem would be first to compute
the inverse of $\textbf{F}^{T}\textbf{F}$, then calculate
$\textbf{F}^{T}\textbf{g}_{l}$ and finally multiply both the entities
together. Notably, the last two operations are matrix-vector products that can
be trivially parallelized. The expensive inverse needs to be done only once.

Now, contrast the above with the following non-linear least squares set-up

$\displaystyle\min_{\textbf{s}_{l}}\|\textbf{f}(\textbf{s}_{l})\|_{2}^{2}\approx\min_{\textbf{s}_{l}}\|\textbf{F}_{l}\textbf{s}_{l}-\textbf{g}_{l}\|_{2}$
(3)

The solution process begins by computing the Jacobian of f around a given
guess solutions $\textbf{s}_{l}$ to obtain a linear least squares
approximation. Importantly, since each $\textbf{s}_{l}$ would be different,
the matrix $\textbf{F}_{l}$ will vary across problem instances. Thus,
computing the solution would first require forming
$\textbf{F}_{l}^{T}\textbf{F}_{l}$ followed by computing inverses (or just
factorization) of each of these. Furthermore, these computations need to be
refined at each iteration of the non-linear least squares.

The above example illustrates the relative difficulty of parallelizing the
per-iteration computations of a non-linear least squares problem compared to
the special linear setting presented in (2). Many trajectory optimizations are
indeed formulated as non-linear least-squares [teb_planner], and thus, they
inherit the same bottlenecks discussed above towards parallelization. We note
that in the non-linear setting, instead of parallelizing per-iteration
operations, it will be more reasonable to solve each non-linear least squares
in full in parallel CPU threads.

The core feature of our batch optimizer is that its most computationally heavy
part has the same structure as the linear least-squares set-up of (2). Thus
our parallelization effort essentially reduces to computing batch matrix-
vector products.

###  III-B Building Blocks

Reformulating the Collision Avoidance Constraints: We rephrase the quadratic
collision avoidance constraints (1e) into the form
$\textbf{f}_{c,l}=\textbf{0}$ based on our prior works [aks_iros21].

$\displaystyle\textbf{f}_{c,l}=\left\\{\begin{array}[]{lcr}x_{l}(t)-\xi_{xj}(t)-ad_{j,l}(t)\cos\alpha_{j,l}(t),\forall
j,t\\\ y_{l}(t)-\xi_{yj}(t)-bd_{j,l}(t)\sin\alpha_{j,l}(t),\forall j,t\\\
\end{array}\right\\}$ (6)

As evident, $\textbf{f}_{c,l}$ resembles a polar representation of the
Euclidean distance between the robot and the obstacle with the variables
$\alpha_{j,l}(t)$ being the angle of the line of sight vector connecting the
ego-vehicle and the $j^{th}$ obstacle. The variable $d_{j,l}(t)$ is the ratio
of the length of the line-of-sight vector to minimum separation distance
required for collision avoidance. Note that these variables will be different
for each problem in the batch and the additional subscript $l$ has been
introduced to represent that fact. Collision avoidance can be enforced by
ensuring $d_{j,l}(t)\geq 1$. It should be noted that in (6),
$\alpha_{j,l}(t),d_{j,l}(t)$ are unknown variables that are obtained by our
optimizer along with other trajectory variables.

Reformulating Acceleration Bounds: Typically, the quadratic acceleration
bounds in (1d) are split into separate affine constraints along each motion
direction. However, we reformulate it in the same manner as collision
avoidance constraints. That is we have constraints of the form
$\textbf{f}_{a,l}=\textbf{0}$, where

$\displaystyle\textbf{f}_{a,l}=\left\\{\begin{array}[]{lcr}\ddot{x}_{l}(t)-d_{a,l}(t)\cos\alpha_{a,l}(t)\\\
\ddot{y}_{l}(t)-d_{a,l}(t)\sin\alpha_{a,l}(t)\\\
\end{array}\right\\},,d_{a,l}(t)\leq a_{max}$ (9)

The variables $\alpha_{a,l}(t),d_{a,l}$ are unknown and will be computed by
our optimizer. On the surface, our representation of acceleration bounds seem
more complicated but as shown later, is essential achieve appropriate
computational structure in our batch optimizer.

Trajectory Parameterization The trajectory variables in each of the $l$
instantiations of the problem can be represented in the following manner.

$\displaystyle\begin{bmatrix}x(t_{1}),\dots,x(t_{n})\end{bmatrix}=\textbf{P}\textbf{c}_{x,l},\begin{bmatrix}\psi(t_{1}),\dots,\psi(t_{n})\end{bmatrix}=\textbf{P}\textbf{c}_{\psi,l}$
(10)

where, P is a matrix formed with time-dependent basis functions (e.g
polynomials) and $\textbf{c}_{x,l},\textbf{c}_{\psi,l}$ are the coefficients
associated with the basis functions. Similar expressions can be written for
$y_{l}(t)$ as well. We can also express the derivatives in terms of
$\dot{\textbf{P}},\ddot{\textbf{P}}$. Our representation creates a low
dimensional parametrization of the trajectories and has its own pros and cons.
On one hand, it reduces the size of the optimization problem but on the other
it also leads to less free parameters for the satisfaction of the constraints.
We have observed that the representation (10) has performed well with our
multi-convex approximation-based approach for trajectory optimization for not
only the current work but also prior efforts like [aks_icra20]. However, it is
unclear how well our low dimensional representation will work with standard
SQP or GD-based approaches.

Matrix Representation: Using the trajectory parametrization presented in (10),
we can put constraints (6) and (9) and the non-holonomic constraints (1b) in
the matrix form of (11a), (11b), (11c) respectively.

$\displaystyle\textbf{F}_{o}\textbf{c}_{x,l}=\boldsymbol{\xi}_{x}+a\textbf{d}_{l}\cos\boldsymbol{\alpha}_{l},\textbf{F}_{o}\textbf{c}_{y,l}=\boldsymbol{\xi}_{y}+b\textbf{d}_{l}\sin\boldsymbol{\alpha}_{l}$
(11a)
$\displaystyle\ddot{\textbf{P}}\textbf{c}_{x,l}=\textbf{d}_{a,l}\cos\boldsymbol{\alpha}_{a},\ddot{\textbf{P}}\textbf{c}_{y,l}=\textbf{d}_{a,l}\sin\boldsymbol{\alpha}_{a},$
(11b)
$\displaystyle\dot{\textbf{P}}\textbf{c}_{x,l}=\textbf{v}_{l}\cos\boldsymbol{\textbf{P}\textbf{c}_{\psi,l}},\dot{\textbf{P}}\textbf{c}_{y,l}=\textbf{v}_{l}\sin\boldsymbol{\textbf{P}\textbf{c}_{\psi,l}}$
(11c)

The matrix $\textbf{F}_{o}$ is obtained by stacking the matrix P from (10) $m$
times (the number of obstacles in the environment). The vector
$\boldsymbol{\xi}_{x},\boldsymbol{\xi}_{y}$ is formed by appropriately
stacking $\xi_{xj}(t),\xi_{yj}(t)$ at different time instants and for all the
obstacles. Similar construction is followed to obtain
$\boldsymbol{\alpha}_{l},\boldsymbol{\alpha}_{a,l},\textbf{d}_{l},\boldsymbol{d}_{a,l},\boldsymbol{\psi}_{l},\textbf{v}_{l}$.

###  III-C Multi-Convex Reformulation

Using previous derivations, we substitute (1a)-(1e) by the following:

$\displaystyle\min_{v_{min}\leq\textbf{v}_{l}\leq\textbf{v}_{max}}\frac{1}{2}\textbf{c}_{x,l}^{T}\textbf{Q}\textbf{c}_{x,l}+\frac{1}{2}\textbf{c}_{y,l}^{T}\textbf{Q}\textbf{c}_{y,l}+\frac{1}{2}\textbf{c}_{\psi,l}^{T}\textbf{Q}\textbf{c}_{\psi,l}$
(12a) $\displaystyle\textbf{A}\begin{bmatrix}\textbf{c}_{x,l}\\\
\textbf{c}_{y,l}\end{bmatrix}=\textbf{b}_{l},\textbf{A}\textbf{c}_{\psi,l}=\textbf{b}_{\psi,l}$
(12b) $\displaystyle\textbf{F}\begin{bmatrix}\textbf{c}_{x,l}\\\
\textbf{c}_{y,l}\\\
\end{bmatrix}=\textbf{g}_{l}(\textbf{c}_{\psi,l},\boldsymbol{\alpha}_{l},\boldsymbol{\alpha}_{a,l},\textbf{d}_{l},\textbf{d}_{a,l},\textbf{v}_{l})$
(12c)

$\displaystyle\textbf{F}=\begin{bmatrix}\begin{bmatrix}\textbf{F}_{o}\\\
\ddot{\textbf{P}}\\\ \dot{\textbf{P}}\end{bmatrix}&\textbf{0}\\\
\textbf{0}&\begin{bmatrix}\textbf{F}_{o}\\\ \ddot{\textbf{P}}\\\
\dot{\textbf{P}}\end{bmatrix}\end{bmatrix},\textbf{g}_{l}=\begin{bmatrix}\boldsymbol{\xi}_{x}+a\textbf{d}_{l}\cos\boldsymbol{\alpha}_{l}\\\
\textbf{d}_{a,l}\cos\boldsymbol{\alpha}_{a,l}\\\
\textbf{v}_{l}\cos\textbf{P}\textbf{c}_{\psi,l}\\\
\boldsymbol{\xi}_{y}+b\textbf{d}_{l}\sin\boldsymbol{\alpha}_{l}\\\
\textbf{d}_{a,l}\sin\boldsymbol{\alpha}_{a,l}\\\
\textbf{v}_{l}\sin\textbf{P}\textbf{c}_{\psi,l}\\\ \end{bmatrix},$ (13)

The cost function is a matrix representation of the sum of squared
acceleration term in (1a). The equality constraints (12b) are the matrix
representation of the boundary constraints (1c). We stack all the non-convex
equality constraints in (12c).

######  Remark 1.

The matrices
$(\textbf{Q},\textbf{F},\textbf{F}_{o},\textbf{P},\dot{\textbf{P}},\ddot{\textbf{P}})$
in the reformulated problem (12a)-(12c) do not depend on the batch index $l$.
In other words, they are the same for all the problem instantiations.

Remark 1 highlights the motivation behind choosing the specific representation
of the collision avoidance (6) and acceleration bounds (9).

###  III-D Solution by Alternating Minimization

We solve (12a)-(12c) by relaxing the non-convex equality constraints (12c) as
$l_{2}$ penalties and augmenting them into the cost function.

$f_{xy}(\textbf{c}_{x,l},\textbf{c}_{y,l},\boldsymbol{\lambda}_{x,l},\boldsymbol{\lambda}_{y,l})+f_{\psi}(\textbf{c}_{\psi,l},\boldsymbol{\lambda}_{\psi,l})+\frac{\rho_{xy}}{2}\left\|\textbf{F}\begin{bmatrix}\textbf{c}_{x,l}\\\
\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}_{l}\right\|_{2}^{2}$ (14)
$\displaystyle
f_{xy}=\frac{1}{2}\textbf{c}_{x,l}^{T}\textbf{Q}\textbf{c}_{x,l}+\frac{1}{2}\textbf{c}_{y,l}^{T}\textbf{Q}\textbf{c}_{y,l}-\langle\boldsymbol{\lambda}_{x,l},\textbf{c}_{x,l}\rangle-\langle\boldsymbol{\lambda}_{y,l},\textbf{c}_{y,l}\rangle$
(15a) $\displaystyle
f_{\psi}=\frac{1}{2}\textbf{c}_{\psi,l}^{T}\textbf{Q}\textbf{c}_{\psi,l}-\langle\boldsymbol{\lambda}_{\psi,l},\textbf{c}_{\psi,l}\rangle$
(15b)

Algorithm 1 Alternating Minimization based Batch Non-Holonomic Trajectory
Optimization

1:while $k\leq maxiter$ do
$\displaystyle{{}^{k+1}}(\textbf{c}_{x,l},\textbf{c}_{y,l})=\arg\min_{\textbf{c}_{x,l},\textbf{c}_{y,l}}f(\textbf{c}_{x},\textbf{c}_{y},{{}^{k}}\boldsymbol{\lambda}_{x,l},{{}^{k}}\boldsymbol{\lambda}_{y,l})+\frac{\rho_{xy}}{2}\left\|\textbf{F}\begin{bmatrix}\textbf{c}_{x,l}\\\
\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}({{}^{k}}\boldsymbol{\alpha}_{l},{{}^{k}}\boldsymbol{\alpha}_{a,l},{{}^{k}}\textbf{c}_{\psi,l},{{}^{k}}\textbf{d}_{l},{{}^{k}}\textbf{d}_{a,l},{{}^{k}}\textbf{v}_{l})\right\|_{2}^{2},\textbf{A}\begin{bmatrix}\textbf{c}_{x,l}\\\
\textbf{c}_{y,l}\end{bmatrix}=\textbf{b}_{l}$ (16)
$\displaystyle{{}^{k+1}}\textbf{c}_{\psi}=\arg\min_{\textbf{A}\textbf{c}_{\psi,l}=\textbf{b}_{\psi,l}}f(\textbf{c}_{\psi})+\frac{\rho_{xy}}{2}\|\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}(\textbf{c}_{\psi,l})\|_{2}^{2}=\arg\min_{\textbf{A}\textbf{c}_{\psi,l}=\textbf{b}_{\psi,l}}f(\textbf{c}_{\psi})+\frac{\rho_{xy}}{2}\left\|\begin{matrix}{{}^{k+1}}\dot{\textbf{x}}_{l}-{{}^{k}}\textbf{v}_{l}\cos\textbf{P}\textbf{c}_{\psi,l}\\\
{{}^{k+1}}\dot{\textbf{y}}_{l}-{{}^{k}}\textbf{v}_{l}\sin\textbf{P}\textbf{c}_{\psi,l}\end{matrix}\right\|_{2}^{2}$
$\displaystyle=\arg\min_{\textbf{A}\textbf{c}_{\psi,l}=\textbf{b}_{\psi,l}}f(\textbf{c}_{\psi})+\frac{\rho_{xy}}{2}\|\arctan
2({{}^{k+1}}\dot{\textbf{y}}_{l},{{}^{k+1}}\dot{\textbf{x}}_{l})-\textbf{P}\textbf{c}_{\psi,l}\|_{2}^{2}$
(17)
$\displaystyle{{}^{k+1}}\textbf{v}_{l}=\arg\min_{v_{min}\leq\textbf{v}_{l}\leq\textbf{v}_{max}}\left\|\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\\\
\end{bmatrix}-\textbf{g}({{}^{k+1}}\textbf{c}_{\psi,l},\textbf{v}_{l})\right\|_{2}^{2}=\arg\min_{v_{min}\leq\textbf{v}\leq\textbf{v}_{max}}\left\|\begin{matrix}{{}^{k+1}}\dot{\textbf{x}}_{l}-\textbf{v}_{l}\cos\boldsymbol{\psi}_{l}\\\
{{}^{k+1}}\dot{\textbf{y}}_{l}-\textbf{v}_{l}\sin\boldsymbol{\psi}_{l}\\\
\end{matrix}\right\|_{2}^{2}$ (18)
$\displaystyle{{}^{k+1}}\boldsymbol{\alpha}_{l}=\arg\min_{\boldsymbol{\alpha}_{l}}\left\|\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}\right\|_{2}^{2}=\arg\min_{\boldsymbol{\alpha}_{l}}\left\|\begin{matrix}{{}^{k+1}}\textbf{x}_{l}-\boldsymbol{\xi}_{x}-a\hskip
2.84544pt{{}^{k}}\textbf{d}_{l}\cos\boldsymbol{\alpha}_{l}\\\
{{}^{k+1}}\textbf{y}_{l}-\boldsymbol{\xi}_{y}-b\hskip
2.84544pt{{}^{k}}\textbf{d}_{l}\sin\boldsymbol{\alpha}_{l}\\\
\end{matrix}\right\|_{2}^{2},{{}^{k+1}}\textbf{d}_{l}=\arg\min_{\textbf{d}_{l}\geq
1}\left\|\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}({{}^{k+1}}\boldsymbol{\alpha}_{l})\right\|_{2}^{2}$
(19)
$\displaystyle{{}^{k+1}}\boldsymbol{\alpha}_{a,l}=\arg\min_{\boldsymbol{\alpha}_{a,l}}\left\|\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}({{}^{k}}\textbf{d}_{a,l},\boldsymbol{\alpha}_{a,l})\right\|_{2}^{2},{{}^{k+1}}\textbf{d}_{a,l}=\arg\min_{\textbf{d}_{a,l}\leq
a_{max}}\left\|\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\end{bmatrix}-\textbf{g}({{}^{k+1}}\boldsymbol{\alpha}_{a,l},\textbf{d}_{a,l})\right\|_{2}^{2}$
(20)

2:end while

Note the introduction of so-called Lagrange multipliers
$\boldsymbol{\lambda}_{x,l}$, $\boldsymbol{\lambda}_{y,l}$, and
$\boldsymbol{\lambda}_{\psi,l}$ that play a crucial role in driving the
residuals of the equality constraints to zero [admm_neural].

######  Remark 2.

The augmented cost function (14) is convex in
$(\textbf{c}_{x,l},\textbf{c}_{y,l})$ for a given
$\textbf{c}_{\psi,l},\boldsymbol{\alpha}_{l},\boldsymbol{\alpha}_{a,l},\textbf{d}_{l},\textbf{d}_{a,l}$.
Similarly, it is convex in $\textbf{d}_{l},\textbf{d}_{a,l}$ for a given
$\textbf{c}_{\psi,l},\boldsymbol{\alpha}_{l},\boldsymbol{\alpha}_{a,l},\textbf{c}_{x,l},\textbf{c}_{y,l}$

######  Remark 3.

For a given $(\textbf{c}_{x,l},\textbf{c}_{y,l})$, the cost function (14) is
non-convex in $\textbf{c}_{\psi,l}$ but can be replaced with a simple convex
surrogate from [aks_icra20].

######  Remark 4.

For a given $(\textbf{c}_{x,l},\textbf{c}_{y,l})$, the optimizations over
variables $(\boldsymbol{\alpha}_{l},\boldsymbol{\alpha}_{a,l})$ have a simple
closed form solution.

######  Remark 5.

The augmented Lagrangian based reformulation of the non-convex constraints
ensures that our batch optimization is always feasible. As a result, it can
handle infeasible (e.g with respect to collision avoidance) trajectory
initialization.

Remarks 2 and 3 are precisely the multi-convex structure foreshadowed in the
earlier sections. Moreover, remarks 2-4 highlight why an AM approach would be
suitable: by decomposing the optimization process over separate blocks of
variable, we can leverage the implicit convex structures present in the
problem. The use of Augmented Lagrangian cost (14) in combination with AM
procedure is known as the split-Bregman technique [admm_neural].

The different steps of AM are presented in (16)-(20), wherein the left
superscript $k$ is used to track the value of the variable over different
iterations. For example, ${{}^{k}}\textbf{c}_{x,l}$ represents the value at
iteration $k$ of this specific variable. At each optimization block, only few
specific variables are optimized while the rest are kept fixed at the values
obtained in the previous iteration or the previous step of the same iteration.

###  III-E Analysis

Step (16) This optimization is a convex equality constrained QP that reduces
to solving the following set of linear equations, wherein
$\boldsymbol{\mu}_{x,l},\boldsymbol{\mu}_{y,l}$ are the dual variable
associated with the equality constraints.

$\displaystyle\overbrace{\begin{bmatrix}\textbf{Q}+\rho_{xy}\textbf{F}^{T}\textbf{F}&\textbf{A}^{T}\\\
\textbf{A}&\textbf{0}\end{bmatrix}}^{\textbf{Q}_{xy}}\begin{bmatrix}\textbf{c}_{x,l}\\\
\textbf{c}_{y,l}\\\ \boldsymbol{\mu}_{x,l}\\\ \boldsymbol{\mu}_{y,l}\\\
\end{bmatrix}=\overbrace{\begin{bmatrix}\rho_{xy}\textbf{F}^{T}{{}^{k}}\textbf{g}_{l}+\begin{bmatrix}{{}^{k}}\boldsymbol{\lambda}_{x,l}\\\
{{}^{k}}\boldsymbol{\lambda}_{y,l}\\\ \end{bmatrix}\\\
\textbf{b}_{l}\end{bmatrix}}^{\textbf{q}_{l}}$ (21)

The set of equations (21) computes the solution for the $l^{th}$ instance of
the problem. However, since the left hand side of (21) does not depend on the
batch index $l$, we can compute the solution of the entire batch in one-shot
through (22).

$\begin{bmatrix}\textbf{c}_{x,1},\textbf{c}_{y,1},\boldsymbol{\mu}_{x,1},\boldsymbol{\mu}_{y,1}\\\
\vdots\\\
\textbf{c}_{x,l},\textbf{c}_{y,l},\boldsymbol{\mu}_{x,l},\boldsymbol{\mu}_{y,l}\end{bmatrix}=\left(\textbf{Q}_{xy}^{-1}\begin{bmatrix}\textbf{q}_{1}|\textbf{q}_{2}|\textbf{q}_{1}\dots\textbf{q}_{l}\end{bmatrix}\right)^{T}$
(22)

The major computation cost of (22) stems from obtaining different
$\textbf{F}{{}^{k}}\textbf{g}_{l}$. But it is straightforward to formulate
this operation as one large matrix-vector product and subsequently parallelize
its computation.

Step (17)  As mentioned earlier, optimization over $(\textbf{c}_{\psi,l})$ is
non-convex due to the presence of the non-holonomic penalty (second term).
However, as shown in the last line of (17), for a given
$({{}^{k+1}}\dot{\textbf{x}}_{l}{{}^{k+1}}\dot{\textbf{y}}_{l})$, the non-
convex term can be replaced with a convex surrogate over
$\textbf{c}_{\psi,l}$, thereby reducing our problem to an equality-constrained
QP. The solution process boils down to solving the following set of linear
equations

$\overbrace{\begin{bmatrix}\textbf{Q}+\rho_{xy}\textbf{P}^{T}\textbf{P}&\textbf{A}^{T}\\\
\textbf{A}&\textbf{0}\end{bmatrix}}^{\textbf{Q}_{\psi}}\begin{bmatrix}\textbf{c}_{\psi,l}\\\
\boldsymbol{\mu}_{\psi,l}\end{bmatrix}=\overbrace{\begin{bmatrix}\rho_{xy}\textbf{P}^{T}\arctan
2({{}^{k+1}}\dot{\textbf{y}}_{l},{{}^{k+1}}\dot{\textbf{x}}_{l})+{{}^{k}}\boldsymbol{\lambda}_{\psi,l}\\\
\textbf{b}_{\psi,l}\end{bmatrix}}^{\textbf{q}_{\psi,l}}$ (23)

Similar to the previous step, the left hand side of (23) do not depend on the
batch index $l$ and thus we can compute the solution for the entire batch in
one-shot

$\begin{bmatrix}\textbf{c}_{\psi,1},\boldsymbol{\mu}_{\psi,1}\\\ \vdots\\\
\textbf{c}_{\psi,l},\boldsymbol{\mu}_{\psi,l}\\\
\end{bmatrix}=\left(\textbf{Q}_{\psi}^{-1}\begin{bmatrix}\textbf{q}_{\psi,1}|\textbf{q}_{\psi,2}|\textbf{q}_{\psi,3}\dots\textbf{q}_{\psi,l}\end{bmatrix}\right)^{T}$
(24)

Step (18): For a given
${{}^{k+1}}\textbf{c}_{x,l},{{}^{k+1}}\textbf{c}_{y,l},{{}^{k+1}}\textbf{c}_{\psi,l}$
or alternately
$({{}^{k+1}}\dot{\textbf{x}}_{l},{{}^{k+1}}\dot{\textbf{y}}_{l},{{}^{k+1}}\textbf{c}_{\psi,l})$,
the velocity $v_{l}(t)$ at different time instants can be treated as
independent of each other. In other words, each element of $\textbf{v}_{l}$ is
decoupled and thus the optimization (18) reduces to $n$ parallel single-
variable quadratic programming problems with a closed form solution. For
$v_{min}>0$, the solution is given by (25), wherein $clip(.)$ performs simple
thresholding to satisfy the velocity bounds.

$\displaystyle{{}^{k+1}}\textbf{v}_{l}=clip(\sqrt{{{}^{k+1}}\dot{\textbf{x}}_{l}^{2}+{{}^{k+1}}\dot{\textbf{y}}_{l}^{2}},v_{min},v_{max})$
(25)

First part of Step (19): For a given
$({{}^{k+1}}{\textbf{x}},{{}^{k+1}}{\textbf{y}})$, each element of
$\boldsymbol{\alpha}_{l}$ can be considered to be decoupled from each other.
Thus, the first optimization in (19) separates into $n$ decoupled problems
with the following closed form solution

$\displaystyle{{}^{k+1}}\boldsymbol{\alpha}_{l}=\arctan
2(a({{}^{k+1}}\textbf{y}_{l}-\boldsymbol{\xi}_{y}),b({{}^{k+1}}\textbf{x}_{l}-\boldsymbol{\xi}_{x}))$
(26)

Second part of Step (19):: Similar to previous step, each element of
$\textbf{d}_{l}$ can be considered independent and thus optimization over
$\textbf{d}_{l}$ reduces to $n$ parallel single-variable QP with simple bound
constraints. We can obtain a closed form solution by first solving the
unconstrained problem and then simply clipping the value to lie between
$[0\hskip 2.84544pt1]$.

Step (20): The two optimizations in this step have the same structure as (19)
and thus have a closed form solution.

######  Remark 6.

The solution of optimization (18)-(20) involves only element-wise operations
without any need of computing matrix factorization or inverse. Thus, computing
a batch solution is trivial.

Multiplier Update: The Lagrange multipliers are updated in the following
manner [admm_neural] which can be trivially done over the entire batch in one-
shot.

$({{}^{k+1}}\boldsymbol{\lambda}_{x,l},{{}^{k+1}}\boldsymbol{\lambda}_{y,l})=({{}^{k}}\boldsymbol{\lambda}_{x,l},{{}^{k}}\boldsymbol{\lambda}_{y,l})-\rho_{xy}\textbf{F}^{T}(\textbf{F}\begin{bmatrix}{{}^{k+1}}\textbf{c}_{x,l}\\\
{{}^{k+1}}\textbf{c}_{y,l}\\\ \end{bmatrix}-{{}^{k+1}}\textbf{g})$ (27)
$\displaystyle{{}^{k+1}}\boldsymbol{\lambda}_{\psi,l}={{}^{k}}\boldsymbol{\lambda}_{\psi,l}-\rho_{xy}\textbf{P}^{T}(\arctan
2({{}^{k+1}}\dot{\textbf{y}},{{}^{k+1}}\dot{\textbf{x}})-\textbf{P}{{}^{k+1}}\textbf{c}_{\psi,l})$
(28)

###  III-F Goal Sampling and Meta-Cost

This section provides a goal-sampling procedure for our batch optimizer and a
meta-cost function to rank the resulting trajectories. We consider two typical
scenarios encountered in highway driving.

Cruise Driving: Our first scenario considers driving forward with velocity as
close as possible to a given $v_{cruise}$. Thus, our meta-cost is defined as
simply.

$\displaystyle\sum_{t}(v(t)-v_{cruise})^{2}.$ (29)

The goal position are spread evenly on different lanes, each at a distance of
$v_{cruise}*t_{f}$, where $t_{f}$ is the planning horizon.

Driving with Maximum Speed close to the Right Lane: In this scenario, the ego-
vehicle is required to drive as close as possible to maximum speed $v_{max}$
while being as close to the right-lane. The meta-cost is defined as the
following wherein $y_{rl}$ is the lateral coordinate of the right-lane and
$w_{1}$ and $w_{2}$ are user-defined constants.

$\displaystyle\sum_{t}w_{1}(v(t)-v_{max})^{2}+w_{2}(y(t)-y_{rl})^{2}.$ (30)

The goals are sampled in the following manner. Around $60\%$ of the goals are
placed on the right-lane at different distances. The remaining goals are
spread across different lanes at a distance of $v_{max}*t_{f}$ from the
current position.

##  IV Validation and Benchmarks

Implementation Details:We implemented our batch optimizer in C++ using Eigen
[eigen]. We used $l=11,a=5.6,b=3.1$ in the simulations. For each driving
scenario discussed in the previous section, we created two variants depending
on whether the neighboring vehicles follow the synthetic Intelligent Driver
Model (IDM) or the pre-recorded trajectories from NGSIM data-set (6 different
scenes) [ngsim]. In the IDM data-set, each neighboring vehicle moves parallel
to center-line and just adapt their cruise forward velocity based on the
distance to the vehicles in front. In the NGSIM data set, the neighboring
vehicles executes the pre-recorded trajectories. It is worth reiterating that
our batch MPC and all the baselines have access to only the instantaneous
position and velocity of the neighboring vehicles and not their true
trajectories over the planning horizon.

Optimizer Convergence: Fig.2(a) shows the typical residual curve obtained with
our batch optimizer in one of the MPC cycles. We show the trend for the best
performing trajectory in the batch. We observed that on an average 100
iterations are enough to obtain residuals in the range of $10^{-3}$ for all
the constraints.

Table II: Meta-Cost Values in different Driving Scenarios (Mean/Min/Max). Lower is better Method | Cruise driving (IDM) | Cruise driving (NGSIM) | High-speed driving (IDM) | High-speed driving (NGSIM)
---|---|---|---|---
Standard MPC | 5.41 / 0.0 / 50.97 | 3.260 / 0.0 / 54.3 | 2141.4 / 1360.8 / 2668.75 | 997.30 / 582.59 / 1256.1
Ours | 0.01 / 0.0 / 0.05 | 0.057 / 0.0 / 0.44 | 238.0 / 135.62 / 425.14 | 236.62 / 142.14 / 574.85
ACADO batch size 11 | 0.08 / 0.0 / 0.66 | 0.114 / 0.0 / 0.62 | 381.96 / 183.01 / 1304.08 | 376.29 / 149.93 / 880.42
ACADO batch size 6 | 0.12 / 0.0 / 1.06 | 0.103 / 0.0 / 0.85 | 643.97 / 129.5 / 1379.8 | 403.46 / 169.69 / 646.85
frenet-frame planner | 0.14 / 0.0 / 1.00 | 0.280 / 0.02 / 0.95 | 563.71 / 194.64 / 1276.27 | 640.48 / 323.77 / 1144.0
Table III: Acceleration Magnitudes Across Different Scenarios(Mean/Min/Max). Lower is better for Cruise driving. For high-speed driving, higher linear acceleration is better. Method | Lin. Acc. Cruise driving | Lin. Acc. High-speed driving | Ang. Acc. Cruise driving | Ang. Acc. High-speed driving
---|---|---|---|---
Standard MPC | 0.93 / 0.00 / 2.63 | 0.72 / 0.00 / 1.56 | 0.02 / 0.00 / 0.07 | 0.01 / 0.00 / 0.03
Ours | 0.11 / 0.00 / 0.28 | 0.99 / 0.00 / 1.77 | 0.02 / 0.00 / 0.07 | 0.43 / 0.00 / 0.16
ACADO batch size 11 | 0.39 / 0.00 . 0.99 | 0.60 / 0.00 / 1.69 | 0.03 / 0.00 / 0.08 | 0.04 / 0.00 / 0.12
ACADO batch size 6 | 0.43 / 0.00 / 1.10 | 0.57 / 0.00 / 1.69 | 0.03 / 0.00 / 0.07 | 0.03 / 0.00 / 0.12
frenet-frame planner | 1.21 / 0.00 / 0.45 | 0.19 / 0.00 / 0.48 | 0.15 / 0.00 / 0.44 | 0.13 / 0.00 / 0.40

Baselines: We benchmark against the following baselines

  * •

Standard MPC: We formulate a single batch MPC wherein the meta cost function
is directly embedded into the trajectory optimizer to compute a single locally
optimal trajectory. We use state-of-the-art optimal control framework ACADO
[acado] as the solver for the standard MPC.

  * •

MPC with Batch ACADO: We construct a batch version of ACADO, which solves
several goal directed MPC over parallel CPU threads. This parallelization does
not require any changes to be made in the matrix-algebra on the underlying SQP
solver in ACADO. Thus, we use this set-up to highlight the computation gain
resulting from our batch solver wherein the per-iteration computation itself
vectorizes across problem instances.

  * •

Frenet Frame Planner: We also compare our batch optimizer with trajectory
sampling approach presented in [fernet_planner] which has been extensively
used in the autonomous driving community and inspired similar related
approaches like [behavior_ad].

###  IV-A Benchmarking

Cruise Scenario: Table II quantify the performance in the cruise-driving
scenario obtained with different methods. The first two columns show the
statistics of the velocity residuals $(v(t)-v_{cruise})^{2}$ observed over the
full run of the MPC. The standard MPC performs the worst with a mean residual
of $5.41$ and $3.26$ on the synthetic (IDM) and NGSIM data-set respectively.
The worst-case performance is around $50$ on both data-sets. The performance
is due to the fact in dense traffic scenarios, the standard MPC trajectories
are unable to find a trajectory around the neighboring slow moving vehicles
(see accompanying video).

Our batch optimizer performs the best with a mean residual of $0.01$ and
$0.05$ on the synthetic (IDM) and NGSIM data-set. In comparison, ACADO with a
batch size of 11 and 6 shows comparable performance in terms of mean values.
However, our batch optimizer achieves 12 times improvement over the worst-case
numbers on the synthetic data set. The Frenet-frame planner with mean
residuals of $0.14$ and $0.28$ performed worse than our’s and batch ACADO. Its
worst-case residual, though, is ten times ours.

High Speed Driving with Right-Lane Preference: This scenario has two competing
terms in the meta cost: maximizing forward velocity and minimizing lateral
distance to the right-lane. Thus, we adopt a slightly different analysis then
before. We first compare the meta cost value across different methods and
subsequently show how those translate to the physical metrics. Table II (last
two columns) summarizes the former results. As before, the standard MPC
performs worst while our batch optimizer achieves the lowest meta cost value
on both synthetic (IDM) and NGSIM data-set. We map these cost values to the
achieved forward velocity and lateral distance residual in Table IV. For
clarity, we present the combined results over the two data sets. We observe
interesting trends here as each approach attempts to minimize the meta-cost
value by trading off velocities and lateral-distance residual in their own
way.

Our batch optimizer achieves highest mean forward velocity of $19.26m/s$. The
performance of ACADO with a batch size of 11 is comparable to ours at
$18.77m/s$ while Frenet-frame planner’s value stood substantially lower at
$16.65m/s$.

Our batch optimizer also maintains a smaller lateral distance to the right-
lane than parallelized ACADO. Our mean distance is $4.43m$ and in comparison,
ACADO with batch size 11 managed a distance residual of $5.31m$ on average.
The respective values achieved with Frenet-frame planner is smaller than batch
ACADO and even ours since it choose a smaller velocity to quickly converge to
the right-lane.

Acceleration Effort: Table III presents the acceleration statistics observed
across different driving scenarios. For ease of exposition, we combined the
data obtained on the synthetic (IDM) and NGSIM data-set and present the
overall mean, minimum and maximum. In the cruise driving scenario, the ideal
acceleration is zero as the ego-vehicle is required to maintain a constant
velocity. As can be seen from Table III (first column), our batch optimizer
comes very close to the ideal performance with mean and maximum linear
acceleration of $0.11m/s^{2}$ and $0.28m/s^{2}$ respectively. Both the mean
and maximum value of the angular acceleration values are very close to zero.
All the other approaches perform worse. Intuitively, a low acceleration value
suggests that our batch optimizer could continuously navigate to free space
less obstructed by neighboring vehicles. We highlight this explicitly in the
accompanying video. All other approaches perform substantially worse. For
example, the standard MPC’s mean linear acceleration is over 9 times higher
than ours. ACADO with batch size of 11 and 6 uses around 4 times higher
acceleration magnitudes. Unsurprisingly, Frenet-frame planner performs worst
since it ignores collision-avoidance constraints during trajectory generation
process. Thus, it routinely enters a situation from where it needs to apply
either emergency braking or execute a sharp turn to avoid collisions.

For high-speed driving, we want the ego vehicle to drive at max speed. Thus,
in contrast to cruise-driving, here, large accelerations magnitude is indeed
necessary for task fulfillment. As shown in Table IV, since our batch
optimizer results in trajectories that gives more preference to maximizing
speed, our mean accelerations are also highest. On the other hand, Frenet-
frame planner achieved the lowest velocity and thus its linear acceleration
values are also the lowest.

Computation Time: We now present the most important result that is crucial in
understanding the previous results in the appropriate context. Table V
compares the computation-time of ours with all the baselines obtained on a
$i7-8750$ processor with 16 GB RAM. Our batch optimizer with mean time of
$0.07s$ is around $6\times$ faster than ACADO with a batch size of 11. In
other words, ACADO needs substantially more computation budget to be even
loosely competitive with our optimizer. Frenet-frame planner’s timing is
comparable to ours but as shown earlier, it performs the worst among all the
multi-modal baselines in terms of meta-cost and acceleration effort. Our
Frenet frame planner implementation used on an average 500 samples. Increasing
this number, could improve the performance but at the expense of higher
computation time.

Fig.2 shows how computation time scales with batch size in our optimizer and
batched multi-threaded ACADO. Our optimizer shows a linear increase which can
be understood in the following manner. In eqn. (22), the matrix
$\textbf{Q}_{xy}$ is independent of the batch size and only the matrix-vector
on the r.h.s of (22) increases quadratically with it. The rate of increase can
be made linear with simple parallelization of matrix multiplication. In
contrast, batch ACADO solves full SQP in parallel CPU cores. Such
parallelization efforts have thread synchronization overhead. Since the total
cores in a laptop is typically 4 to 6, at batch size of 11 (or even 6), all
SQP instantiations compete with each other for computing resources.
Unfortunately, that is the best that we can achieve with off-the-shelf
solvers. Since SQP relies on linearization, it isn’t easy to parallelize its
computation at each iteration (recall Section III-A).

Table IV: Metrics for High-speed driving (Mean/Min/Max) Method | Lat. dist. from right-lane | Velocity
---|---|---
Standard MPC | 2.72 / 0.0 / 6.36 | 11.16 / 7.66 / 16.45
Ours | 4.43 / 0.0 / 18.0 | 19.28 / 17.14 / 20.68
ACADO batch size 11 | 5.31 / 0.0 / 18.0 | 18.77 / 15.48 / 20.17
ACADO batch size 6 | 3.5 / 0.0 / 15.43 | 17.64 / 14.10 / 20.56
Frenet-frame planner | 2.42 / 0.0 / 9.34 | 16.65 / 14.0 / 19.74
Table V: MPC computation time[s] (Mean/Min/Max) Method | Cruise driving | High-speed driving
---|---|---
Standard MPC | 0.36 / 0.05 / 0.62 | 0.28 / 0.05 / 0.61
Ours | 0.07 / 0.06 / 0.08 | 0.07 / 0.06 / 0.07
ACADO batch size 11 | 0.36 / 0.10 / 0.61 | 0.44 / 0.20 / 0.67
ACADO batch size 6 (infeasible) | 0.25 / 0.11 / 0.39 | 0.32 / 0.15 / 0.47
Frenet-frame planner | 0.082 / 0.068 / 0.097 | 0.08 / 0.07 / 0.096
Figure 2: (a): Residual trend observed for best performing trajectory in a
batch in one of the MPC cycles. We show the different components of constraint
(12c). (b): Computation time vs Batch Size comparison between our batch
optimizer and batched ACADO [acado] parallelized on multiple CPU threads.

##  V Discussions and Future Work

This paper showed how a standard MPC based on local optimization techniques
cannot generate sophisticated driving behaviors. Past works such as
[fernet_planner] have attempted to capture the multi-modality of autonomous
driving by generating many candidate trajectories to different goals and
ranking them based on a cost-function. Importantly, collision-avoidance and
kinematic constraints were ignored during the trajectory generation process,
leading to poor fulfilment of the given driving tasks. As a potential
solution, we presented a trajectory optimizer that can generate a batch of
solutions in parallel while incorporating all the necessary constraints. We
showed that potential competing baselines based on state-of-the-art optimal
control solver ACADO need up to 6x more computation time budget to produce
comparable solution as ours.

Our meta-costs can capture some of the rule-set described in [rule_set_knoll].
For example, the velocity residual (29) can be easily modified to induce
behaviors where ego-vehicle overtakes only when it leads to some minimum gain
in the forward velocity. Similarly, "keep-right" rule from [rule_set_knoll] is
already incorporated in the meta-cost (30). We present a more detailed
analysis in our supplementary material [github_repo_driving].

Limitations and Possible Workarounds: Our batch optimization and MPC is
capable of handling curved roads and residential driving scenarios (see
accompanying video). However, it might struggle in highly cluttered and
unstructured environments like parking lot that might require a large heading
change. A possible workaround is to model the ego-vehicle geometry as a
combination of circles, thus explicitly bringing the heading angle in the
collision avoidance model [aks_icra20]. Our preliminary results have shown
that multi-circle approximation preserves the core batch structure at least
for the holonomic robots [aks_gpu_batch_2022].

Our batch optimizer structure, specifically the update rules (22),(24) is
expected to be preserved for the more complex bi-cycle model of the ego-
vehicle and we present a sketch of that derivation in the supplementary
material [github_repo_driving]. However, extension to include dynamics and
tire forces will require major overhaul of the matrix algebra.

Future Work: Our batch optimizer can be an attractive policy class for
reinforcement learning algorithms. It can ensure safe exploration and thus
would be beneficial for deep-Q learning based approaches. Our batch MPC can
also be used to generate supervision data for imitation learning algorithms
like [chauffeurnet]. Our future efforts are geared towards these directions.


2021

These authors contributed equally to this work.

These authors contributed equally to this work.

1]Amrita Center for Cybersecurity Systems & Networks, Amrita Vishwa
Vidyapeetham, Amritapuri, India

2]Department of Gastroenterology, Amrita Institute of Medical Sciences, Kochi,
India

# Automatic Endoscopic Ultrasound Station Recognition with Limited Data

Abhijit Ramesh  K N Anantha Nandanan  Nikhil Boggavarapu  Priya Nair  Gilad
Gressel<EMAIL_ADDRESS>[  [

###### Abstract

Pancreatic cancer is a lethal form of cancer that significantly contributes to
cancer-related deaths worldwide. Early detection is essential to improve
patient prognosis and survival rates. Despite advances in medical imaging
techniques, pancreatic cancer remains a challenging disease to detect.
Endoscopic ultrasound (EUS) is the most effective diagnostic tool for
detecting pancreatic cancer. However, it requires expert interpretation of
complex ultrasound images to complete a reliable patient scan. To obtain
complete imaging of the pancreas, practitioners must learn to guide the
endoscope into multiple “EUS stations” (anatomical locations), which provide
different views of the pancreas. This is a difficult skill to learn, involving
over 225 proctored procedures with the support of an experienced doctor. We
build an AI-assisted tool that utilizes deep learning techniques to identify
these stations of the stomach in real time during EUS procedures. This
computer-assisted diagnostic (CAD) will help train doctors more efficiently.
Historically, the challenge faced in developing such a tool has been the
amount of retrospective labeling required by trained clinicians. To solve
this, we developed an open-source user-friendly labeling web app that
streamlines the process of annotating stations _during_ the EUS procedure with
minimal effort from the clinicians. Our research shows that employing only 43
procedures with no hyperparameter fine-tuning obtained a balanced accuracy of
90%, comparable to the current state of the art. In addition, we employ Grad-
CAM, a visualization technology that provides clinicians with interpretable
and explainable visualizations.

###### keywords:

station classification, Endoscopic ultrasound (EUS), convolutional neural
network, pancreatic cancer

##  1 Introduction

Pancreatic cancer is the seventh leading cause of cancer-related deaths
worldwide [3]. Early detection is essential to improve the prognosis and
survival rate of patients with pancreatic cancer. Despite advances in imaging
technology, the survival rate remains low, with a reported 12% survival rate
[19]. This is because pancreatic cancer is often asymptomatic until it reaches
an advanced stage, making early detection unlikely.

Magnetic Resonance Imaging (MRI), Computed Tomography (CT) scans, Endoscopic
Ultrasound (EUS), and Positron Emission Tomography (PET) scans are different
medical imaging techniques used for diagnosing pancreatic cancer. Among these
imaging techniques, EUS uses high-frequency ultrasound waves to produce
detailed images of the pancreas and surrounding organs. EUS is considered the
most effective method for detecting early pancreatic cancer because it
provides the most accurate visual of the size, location, and extent of the
tumour [7]. A significant advantage of EUS is its capacity to detect very
small tumors, as small as 2mm - 3mm in size; in comparison, CT and MRI can
only detect tumors larger than 1cm. The EUS procedure demands a high level of
expertise and experience, involving over 225 proctored procedures to be done
before being assessed for competency [5]. In other forms of US imaging, the
probe location is fixed and easy to control because it is in the examiners
hand. However, during EUS, clinicians trained in endoscopy need to interpret
complex ultrasound images in real-time using a flexible probe in a constantly
moving environment. To obtain complete imaging of the pancreas, practitioners
must learn to guide the endoscope into multiple “EUS stations” (anatomical
locations), which provide different views of the pancreas. The recognition of
the EUS stations is crucial to the EUS procedure as it enables targeted
biopsies, accurate diagnosis, and aids in further surgical planning and
monitoring.

In order to assist doctors in learning the EUS procedure, previous studies
have demonstrated the feasibility of computer-aided diagnostic (CAD) systems
that use deep learning techniques in order to identify the pancreatic stations
and whether or not the tumor is cancerous [6, 16]. However, these studies
required retrospective annotated data from expert clinicians, increasing the
clinician’s workload. To solve this, we have developed and open-sourced a
labeling application that streamlines the process for the Endoscopist to
annotate the pancreas station during the EUS procedure, adding nearly no
effort to the clinician. This type of “real-time” labeling is successful
because we train our CAD system on all the data found in the video. We do not
require the Endoscopist to only select gold-standard images and manually
remove difficult images.

We build an explainable AI-assisted tool to help Endoscopists identify
different stations during the procedure. Incorporating this AI-assisted tool
can improve the accuracy of diagnoses and decision-making in treating
pancreatic cancer. Importantly, we demonstrate that a state-of-the-art system
can be built with limited data and little labeling effort from the clinicians.
By reducing the data requirements for training models, we aim to democratize
these complex CAD systems. Our study utilized only 43 EUS procedures,
accounting for approximately 10-15% of the data utilized in other related
studies [24, 13]. We leverage preprocessing techniques on the EUS images to
enhance image quality, thereby improving the overall performance of the model
[10]. We also incorporate Grad-CAM visualizations to provide insight into the
decision-making process of deep learning models, thus producing an explainable
CAD system [20, 17]. This also enables the CAD system to be used as an offline
training program, where doctors can practice identifying stations
retrospectively. Our experiments show that using preprocessing techniques, we
achieve an accuracy of 90%, and without using the preprocessing techniques, we
achieve an accuracy of 87.6%.

We make the following novel contributions:

  * •

Achieve a balanced accuracy of 90% with only preprocessing the input and zero
fine-tuning on a small dataset; this is comparable to state-of-the-art that
uses transfer learning, fine-tuning, and larger datasets [23, 24, 13].

  * •

Utilize Grad-CAM to provide explainability to the physicians during
procedures.

  * •

Develop and open-source an EUS labeling application that allows clinicians to
annotate station and FNA timestamps during pancreas ultrasound procedures.
Github Repository.

##  2 Related Work

In recent years, deep learning algorithms have shown great promise in
detecting various diseases from medical imaging data, including ultrasound
images [23]. Studies have shown that deep learning models perform comparably
to or better than human healthcare professionals (HCPs) in most cases,
particularly in the detection of diseases such as skin cancer, breast cancer,
lung cancer, and diabetic retinopathy. These findings suggest that deep
learning algorithms may be effective in the task of pancreas station
classification using ultrasound images [12].

Yao et al. [23], proposed a framework consisting of a bile duct(BD)
segmentation network and station recognition network to classify biliary tract
diseases in EUS images. The framework achieved a classification accuracy of
93.3% for BD stations and an F1 score of 0.77 for segmentation for the
internal validation set. For classification, the model attained an accuracy of
82.6% for the external validation set. This highlights the potential for deep
learning in the context of pancreatic cancer diagnosis and staging.

Several other studies have investigated deep-learning techniques for
pancreatic station classification. Zhang et al. [24] proposed a transfer
learning approach using a pre-trained ResNet model to classify six stations of
the pancreas. The authors achieved an accuracy of 82.4% on the external
dataset and an accuracy of 86.2% on the per-frame video test. The analysis
involved selecting the best images from 311 videos. In our study with 43
videos, less than 15% of their dataset, we achieved an accuracy of 90% on the
test dataset.

LU et al. [13] proposed a two-stage framework where a convolutional neural
network (CNN) was first used to detect the pancreas in the EUS images. Then, a
region-based CNN was trained to classify the pancreas into different EUS
stations. The authors achieved an accuracy of 95.6% in station recognition. To
conduct their investigation, the authors utilized a dataset consisting of 269
procedures, roughly six times our dataset.

Building on this, Jarmillo et al. [9] proposed a novel approach for
automatically detecting pancreatic tumors in endoscopic ultrasound (EUS)
videos using transfer learning. The authors used pre-trained CNN models to
classify cancerous pancreatic tumors vs. benign pancreatic tumors and achieved
an accuracy of 93.2. The CNN was trained on a dataset of 66,249 EUS images
from 55 videos and was evaluated on a test set of 22,097 images from 18
videos. They used pre-processing techniques to remove noise and enhance the
tumor region, resulting in improved accuracy and reduced variability in image
quality.

None of the prior studies systematically experimented with different
preprocessing techniques in order to enhance the performance of deep learning
models on EUS datasets. All previous work relied on large datasets
retrospectively annotated by human clinicians, making scalability difficult.

Despite working with a smaller dataset, comprising just approximately 10-15%
of the videos used in comparative studies, we achieved a balanced accuracy of
90% for the test set. Our data labeling process did not involve manually
annotating “gold standard” images, making it less time-consuming compared to
previous approaches. We also introduce the Grad-CAM visualization technique to
provide clinicians with transparent and explainable results. This allows for a
better understanding of the model’s decision-making process and facilitates
trust for medical professionals.

Table 1: Comparison of EUS Datasets for Station Classification

Paper | N. Patients | N. Images | Performance | Hyperparameter Tuning
---|---|---|---|---
[13] | 269 | 18,061 | 94.1% | Not Specified
[24] | 311 | 19,486 | 94.2% | Fine-tuning
[9] | 55 | 66,249 | 93.2% | Grid search
[6] | 41 | 179,092 | 66.8% | Not Specified
Our Work | 43 | 16,081 | 90% | None

###  2.1 Preprocessing methods

(a) CLAHE

(b) Quantile Capping

(c) Gaussian Smoothing

(d) Fourier Transform

(e) Denoising

Figure 1: Transformed EUS Image under different Preprocessing techniques

As with all machine learning, enhancing data quality will result in improved
performance. Due to the nature of video recording, artifacts, and additional
noise are often captured in a recording, which degrades the quality of the
video. We employ image preprocessing techniques to increase contrast, remove
noise, and smoothen/blur the image.

  * •

Contrast-Limited Adaptive Histogram Equalization: Contrast-limited adaptive
histogram equalization (CLAHE) is an image enhancement technique that is
widely used to improve the contrast and brightness of images [14]. It is a
variant of adaptive histogram equalization (AHE), a non-linear contrast
enhancement method. The idea behind histogram equalization is to transform an
image’s intensity histogram such that the distribution of the pixel
intensities is spread out more evenly across the intensity range, thereby
enhancing the contrast of the image. However, histogram equalization may
produce undesirable results, such as noise amplification in regions of low
contrast. CLAHE was introduced to address this limitation by dividing the
image into smaller regions and equalizing the histogram separately for each
region. This results in a more localized contrast enhancement, which can help
preserve the details of the image [18].

  * •

Gaussian Smoothing: Image filtering techniques, such as Gaussian smoothing,
are commonly used in image denoising to reduce high-frequency noise while
preserving image edges and structures. Gaussian smoothing involves convolving
an image with a Gaussian kernel, a bell-shaped function that assigns weights
to neighbour pixels based on their distance from the center pixel. The
smoothing effect of the Gaussian kernel is determined by its standard
deviation, with a higher standard deviation resulting in a more significant
blur. Xiao et al. [22] presented a study on developing an automatic brain MRI
segmentation scheme with Gaussian smoothing. The paper highlighted the
significance of Gaussian smoothing in medical imaging, showcasing its ability
to preprocess images, effectively reducing noise, and enhancing the clarity of
the features in the image.

  * •

Quantile Capping: Quantile capping is a technique used in image processing to
limit the dynamic range of an image by capping the extreme pixel values. This
is done by finding the upper and lower quantiles of the pixel intensity
distribution and capping the values outside this range. This process can help
improve the visual quality of an image and enhance its contrast. This
technique converts the features into a normal distribution, spreading out the
most common values for a particular feature. It also lowers the impact of
marginal outliers.

  * •

Denoising: In image processing, denoising techniques are used to remove
unwanted noise from images. Non-local Means-Based Denoising(NLM) algorithm is
one such image denoising algorithm that utilizes the self-similarity of images
to remove noises, while still preserving the important image artifacts. it
works by comparing each pixel in the image to all other pixels, and then the
average of these similar pixels is used to get the effective denoised pixel
value. Heo et al. [8] undertook a systematic review to determine the
effectiveness on using NLM algorithm for denoising (Magnetic Resonance)MR
images. The study shows that not only was it effective at removing noises from
the MR images while still keeping the important image artifacts, but also
outperformed other denoising algorithms in terms of peak signal-to-noise
ratio(PSNR) value, which is quality meansurement in which higher the PSNR
value, better the image quality.

  * •

Fourier Transform: The Fourier Transform is a mathematical tool for
decomposing a signal into its frequency components. In image processing, the
Fourier Transform converts an image from its pixel-based spatial domain to its
frequency-based domain. This transformation helps analyze the image in terms
of its frequency content. Working in the frequency domain has several
advantages, such as filtering or enhancing specific frequencies, which can be
performed more efficiently. Fourier Transform is used in image compression and
for the removal of noise in an image. To remove noise in an image, a mask is
applied to the transformed image to suppress the noise frequency components
while preserving the desired image information. Various strategies can be
employed to design effective masks for noise removal in images. This includes
utilizing masks of different shapes and sizes, targeting specific regions of
the transformed image, or applying thresholding techniques to identify and
suppress noise-containing frequencies. After applying a filter to remove noise
frequency components, the inverse Fourier Transform transforms the modified
frequency-based image into the spatial domain.

###  2.2 CNN

The convolutional neural network (CNN) [1] is widely used and considered the
most effective for medical image classification tasks [2]. The CNN is a
powerful feature extractor; therefore, it can be used to classify medical
images and avoid complex and expensive feature engineering.

We used three classical convolutional neural network architectures: ResNet,
DenseNet, and EfficientNet. A description of these architectures is presented
below.

  * •

ResNet: This architecture solves a common problem in deep learning called the
vanishing gradient problem. This problem occurs when deep neural networks have
many layers, which makes it difficult for them to learn effectively. ResNets
addresses this problem by using a technique called skip connections. The skip
connection connects non-contiguous layers using a direct connection. These
connections act like shortcuts, allowing information to flow easily through
the network. By adding these shortcuts, the network can train deeper and
increase performance on classification tasks.

The number of trainable parameters of ResNet18 and ResNet34 is 11M and 63.5M,
respectively.

  * •

DenseNet: This architecture uses dense connections between layers through
Dense Blocks, which directly connect all layers (with matching feature-map
sizes) to one another. Each layer obtains additional inputs from all preceding
layers and passes its feature maps on to all subsequent layers. These dense
connections help the model gain collective knowledge with fewer layers (with
fewer parameters), which helps the model learn faster.

The number of trainable parameters of DenseNet121, DenseNet161 and DenseNet201
is 8M, 29M, and 20M, respectively.

  * •

EfficientNet: This architecture manages the tradeoff between accuracy and
computational cost. It achieves this by using a compound scaling technique
that scales the network’s depth, width, and resolution in a principled manner.
Unlike conventional practice that uses arbitrary values to scale these
factors. The EfficientNet scaling method uniformly scales network width,
depth, and resolution with a set of fixed scaling coefficients $\alpha$,
$\beta$, and $\gamma$. These scaling coefficients are determined by a small
grid search on the original small model [21]. The trainable parameters of
EfficientNetB0 and EfficientNetB3 are 5.3M and 12M, respectively.

##  3 Methodology

###  3.1 Dataset

(a) Label Tool Screen 1

(b) Label Tool Screen 2

(c) Label Tool Screen 3

Figure 2: Screenshots of the Label Tool

Endoscopists from Amrita Hospitals used our open-source label tool shown in
Fig 2 to annotate the timestamps during the EUS procedure. This eliminates the
need for retrospective labelling. The station timestamps from the app
correlate with the endoscopy machine’s screen capture, allowing us to extract
the corresponding frames for a station.

We used 43 clinical videos, which were captured at 24 frames per second (FPS),
and their timestamps of Station 1, Station 2, and Station 3 of the Endoscopy
procedure collected over three months as our dataset for this study. This
dataset is cleaned and prepared for training using our data preparation step.
3.2. The images were extracted from EUS videos at a rate of 1FPS. Previously,
we experimented with different frame rates ranging from 1 to 24, but we
settled on 1FPS because other options reduced the model’s performance. The
dataset comprises 15,545 images, divided as 2,242 for testing and 13303 for
training across all three stations, as shown in Table 2.

We didn’t manually choose or use a special set of ’gold standard’ images for
training our model. The similarity between nearby interval frames means that
testing may involve predicting the class of a patient already in the training
set, making accurate assessment difficult. To address this, patient images
were not mixed across splits. Instead, the dataset was split based on patients
while maintaining balanced proportions of station images in the train, test,
and validation sets. This approach evaluates the model’s generalization
capabilities on unseen images from different patients.

Table 2: Summary of EUS Dataset Station | Train | Test
---|---|---
Station 1 | 4,179 | 744
Station 2 | 5,602 | 830
Station 3 | 3,522 | 668
Total | 13,303 | 2,242

(a) Station 1

(b) Station 2

(c) Station 3

Figure 3: Images of different stations in EUS

###  3.2 Dataset Preparation

We utilized the EUS videos of patients to generate the dataset as mentioned in
the Section 3.1. We identified four distinct types of noises represented in
Fig 4 \- GUI images, white light, green pointers (used by doctors to point on
the screen in real-time), and blackened images.

(a) GUI image

(b) white light image

(c) Green pointer image

(d) Blackened image

Figure 4: Different types of noises in EUS videos

We performed a set of data-cleaning steps on the extracted frames. First, we
removed the GUI and pink images using two histogram comparison techniques -
histogram comparison intersection [11] and histogram comparison Bhattacharya
[4]. Histogram comparison intersection measures the overlap between two
histograms, while histogram comparison Bhattacharya measures the distance
between them. These techniques are chosen because they provide a comprehensive
understanding of the similarity between two histograms of two images. It is
worth noting that the specific threshold values we have used in our processes
were determined through trial and error on our dataset. Thus, they may not be
standard and may need to be adjusted for different datasets. We compared the
extracted images to a reference image. If the histogram comparison
intersection value is less than or equal to 1.031 and the histogram comparison
Bhattacharya value is more than or equal to 0.95, we detected the pink photos.
If the histogram comparison intersection value is larger than or equal to 1.42
and the histogram comparison Bhattacharya value is less than or equal to 0.18,
we remove the GUI image.

Next, we removed the blackened images by computing the average intensity
values of all pixels in the images. If the average pixel value is less than
the threshold of 12, then that image is deemed blackened and removed. Finally,
the images containing green pointers were replaced using an image processing
technique called inpainting [15].

To further enhance the quality and contrast of the extracted frames, we
experimented with several image enhancement algorithms, including Contrast-
limited adaptive histogram equalization(CLAHE), Gaussian smoothing, denoising,
quantile capping, and Fourier Transform as discussed in section 2.1.

The images were normalized and standardized using the mean and standard
deviation in the training set. This is to prevent any data leakage and
promotes robust generalization to unseen data. The architecture diagram of the
proposed framework system is displayed in Fig 5.

Figure 5: Diagram of proposed framework system

###  3.3 Performance Measures

We use balanced accuracy as the primary evaluation metric and weighted
precision and weighted recall as the secondary evaluation metrics. We choose
balanced accuracy as all mistakes are equally weighted; that is, all mistakes
are equally important. Balanced accuracy considers the model’s accuracy in
each class while also considering the number of images in each class. This
ensures that the performance measure is not skewed by the disproportionate
number of images in each class.

$\text{balanced
accuracy}=\frac{1}{k}\sum_{i=1}^{k}\frac{TP_{i}}{TP_{i}+FN_{i}}$ (1)

where,

$TP_{i}$: number of true positives for class $i$

$FN_{i}$: number of false negatives for class $i$

$k$: total number of classes

The use of weighted precision and weighted recall allows us to account for the
imbalance in the dataset by giving more weight to the classes with more
images. This enables us to accurately measure the precision and recall for
each class, taking into consideration the proportion of images in each class.

$precision_{weighted}=\frac{n_{1}P_{1}+n_{2}P_{2}+n_{3}P_{3}}{n_{1}+n_{2}+n_{3}}$
(2)

$recall_{weighted}=\frac{n_{1}R_{1}+n_{2}R_{2}+n_{3}R_{3}}{n_{1}+n_{2}+n_{3}}$
(3)

where,

$P_{i}$: precision for class $i$

$R_{i}$: recall for class $i$

$n_{i}$: number of instances of class $i$

##  4 Results and Discussion

Table 3: Performance Comparisons, best Balanced Accuracy(BA) is bolded

| Resnet18 | Efficientnet_b0 | Densenet161
---|---|---|---
Preprocessing | BA | Precision | Recall | BA | Precision | Recall | BA | Precision | Recall
NO-PRE | 85.6 | 85.7 | 85.6 | 82.2 | 82.3 | 82.1 | 87.6 | 88.5 | 87.7
CLAHE | 84.9 | 85.8 | 85.3 | 80.3 | 80.8 | 80.5 | 80.7 | 82.3 | 81.0
DENOISING | 83.9 | 84.2 | 84.1 | 80.8 | 81.4 | 81.0 | 88.7 | 90.0 | 88.9
GAUSSIAN Smoothing | 84.9 | 86.3 | 85.0 | 81.1 | 81.3 | 81.2 | 90.0 | 90.1 | 90.0
QUANTILE CAP | 75.9 | 76.5 | 76.0 | 77.3 | 78.0 | 77.5 | 59.8 | 70.0 | 61.6
FFT-Normal | 84.4 | 85.8 | 84.6 | 80.2 | 80.4 | 80.2 | 88.5 | 90.0 | 88.8

The results in Table 3 summarise the outcomes of our experiments. The results
of ResNet18, DenseNet161, and EfficientNetB0 trained on the datasets of
preprocessing techniques are shown. DenseNet161, trained on a dataset
preprocessed with Gaussian Smoothing, achieved the highest performance,
attaining a balanced accuracy of 90%. One surprising result is that models
without preprocessing perform only 2.5% lower than the best-performing model.
This has a significant downstream implication, as it provides simplicity in
the workflow and deployment of the model. Without any complex preprocessing
steps, it will increase the inference time of the CAD system during live EUS
procedures.

Furthermore, the preprocessing techniques Gaussian smoothing, denoising, and
FFT outperformed the baseline by 90%, 88.7%, and 88.53%, respectively.
Quantile capping and CLAHE, on the other hand, performed poorly on the
dataset, with balanced accuracy of 77.39% and 84.87%, respectively. This
outcome emphasizes the importance of evaluating the applicability of
preprocessing techniques on a case-by-case basis in medical imaging.

###  4.1 Qualitative Analysis

We conducted a qualitative analysis to evaluate the effectiveness of our deep-
learning models in differentiating stations in the EUS procedure. We utilized
the Grad-CAM technique to visualize the regions of interest (ROIs) that our
model used to make predictions. Grad-CAM is a technique used in deep learning
that produces visual explanations of the decision-making process of a
convolutional neural network by highlighting the regions of an input image
that are most important for the network’s predictions.

According to our findings, the features emphasized by our deep learning model
to classify the EUS images are consistent with the reference points used by
our expert doctor. In essence, the model’s attention mechanism appears to be
focused on the same visual signals that expert doctors use in EUS recordings
to diagnose anomalies.

(a) Station-1

(b) Station-2

(c) Station-3

Figure 6: Grad-CAM visualisation on EUS

##  5 Conclusion

In this work, we demonstrate that it is possible to build a CAD system to help
train doctors in EUS station identification with frugal resources. We created
an open-source labeling tool to annotate the timestamps of pancreas stations
during the EUS procedure. This form of annotation happens in real-time during
the procedure, requiring little extra effort on the behalf of the clinicians.
We then demonstrate that it is possible to achieve state-of-the-art results
with limited data (15% of other studies). Notably the dataset does not have
any hand annotations by doctors for gold-standard images. Instead, we train
our models on the raw video footage with only the time stamps for station
identification. We tested three models with a number of preprocessing
techniques and found the best result was DenseNet161 with an accuracy of 90%.
Furthermore, the result also shows that without preprocessing, the model
achieved an accuracy of 87.6%. It is desirable to drop preprocessing as it
will greatly increase the speed of inference for the real-time application of
the CAD.

Our results thus provide proof that a simplified approach to obtaining
annotated EUS videos and using a small dataset, combined with basic deep
learning techniques, can yield competitive performance. We believe in time and
collaboration, we can reach much higher capacity datasets and obtain
significantly improved models.

##  6 Declarations

###  6.1 Data Availability

The datasets generated during and/or analyzed during the current study are not
publicly available as it is sensitive data for which we do not have consent
for public distribution. Those wishing to collaborate formally are requested
to email.

###  6.2 Conflict of Interest

The authors have no relevant financial or nonfinancial interests to disclose.

## References

  *   * Aloysius and Geetha [2017] Aloysius N, Geetha M (2017) A review on deep convolutional neural networks. In: 2017 international conference on communication and signal processing (ICCSP), IEEE, pp 0588–0592
  * Anwar et al [2018] Anwar SM, Majid M, Qayyum A, et al (2018) Medical image analysis using convolutional neural networks: a review. Journal of medical systems 42:1–13
  * Bray et al [2018] Bray F, Ferlay J, Soerjomataram I, et al (2018) Global cancer statistics 2018: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 68(6):394–424
  * Dalal and Triggs [2005] Dalal N, Triggs B (2005) Histograms of oriented gradients for human detection. In: 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), Ieee, pp 886–893
  * Faulx et al [2017] Faulx AL, Lightdale JR, Acosta RD, et al (2017) Guidelines for privileging, credentialing, and proctoring to perform gi endoscopy. Gastrointestinal endoscopy 85(2):273–281
  * Fleurentin et al [2022] Fleurentin A, Mazellier JP, Meyer A, et al (2022) Automatic pancreas anatomical part detection in endoscopic ultrasound videos. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization pp 1–7
  * Gonzalo-Marin et al [2014] Gonzalo-Marin J, Vila JJ, Perez-Miranda M (2014) Role of endoscopic ultrasound in the diagnosis of pancreatic cancer. World journal of gastrointestinal oncology 6(9):360
  * Heo et al [2020] Heo YC, Kim K, Lee Y (2020) Image denoising using non-local means (nlm) approach in magnetic resonance (mr) imaging: a systematic review. Applied Sciences 10(20):7028
  * Jaramillo et al [2022] Jaramillo M, Ruano J, Gómez M, et al (2022) Automatic detection of pancreatic tumors in endoscopic ultrasound videos using deep learning techniques. In: Medical Imaging 2022: Ultrasonic Imaging and Tomography, SPIE, pp 106–115
  * Keerthi and Santhi [2023] Keerthi S, Santhi P (2023) Precise multi-class classification of brain tumor via optimization based relevance vector machine. Intelligent Automation & Soft Computing 36(1)
  * Lee et al [2005] Lee S, Xin JH, Westland S (2005) Evaluation of image similarity by histogram intersection. Color Research & Application: Endorsed by Inter-Society Color Council, The Colour Group (Great Britain), Canadian Society for Color, Color Science Association of Japan, Dutch Society for the Study of Color, The Swedish Colour Centre Foundation, Colour Society of Australia, Centre Français de la Couleur 30(4):265–274
  * Liu et al [2019] Liu X, Faes L, Kale AU, et al (2019) A comparison of deep learning performance against health-care professionals in detecting diseases from medical imaging: a systematic review and meta-analysis. The lancet digital health 1(6):e271–e297
  * LU et al [2021] LU Z, WU H, YAO L, et al (2021) A station recognition and pancreatic segmentation system in endoscopic ultrasonography based on deep learning. Chinese Journal of Digestive Endoscopy pp 778–782
  * Mohan et al [2020] Mohan M, Nair LS, et al (2020) Fuzzy c-means segmentation on enhanced mammograms using clahe and fourth order complex diffusion. In: 2020 Fourth International Conference on Computing Methodologies and Communication (ICCMC), IEEE, pp 647–651
  * OpenCV [Year] OpenCV (Year) Inpainting tutorial - opencv documentation. Online, URL https://docs.opencv.org/3.4/df/d3d/tutorial_py_inpainting.html
  * Osareh and Shadgar [2011] Osareh A, Shadgar B (2011) A computer aided diagnosis system for breast cancer. International Journal of Computer Science Issues (IJCSI) 8(2):233
  * Panwar et al [2020] Panwar H, Gupta P, Siddiqui MK, et al (2020) A deep learning and grad-cam based color visualization approach for fast detection of covid-19 cases using chest x-ray and ct-scan images. Chaos, Solitons & Fractals 140:110,190
  * Pizer [1990] Pizer SM (1990) Contrast-limited adaptive histogram equalization: Speed and effectiveness stephen m. pizer, r. eugene johnston, james p. ericksen, bonnie c. yankaskas, keith e. muller medical image display research group. In: Proceedings of the first conference on visualization in biomedical computing, Atlanta, Georgia, p 1
  * Rawla et al [2019] Rawla P, Sunkara T, Gaduputi V (2019) Epidemiology of pancreatic cancer: global trends, etiology and risk factors. World journal of oncology 10(1):10–27
  * Selvaraju et al [2017] Selvaraju RR, Cogswell M, Das A, et al (2017) Grad-cam: Visual explanations from deep networks via gradient-based localization. In: Proceedings of the IEEE international conference on computer vision, pp 618–626
  * Tan and Le [2019] Tan M, Le Q (2019) Efficientnet: Rethinking model scaling for convolutional neural networks. In: International conference on machine learning, PMLR, pp 6105–6114
  * Xiao et al [2010] Xiao K, Ho SH, Bargiela A (2010) Automatic brain mri segmentation scheme based on feature weighting factors selection on fuzzy c-means clustering algorithms with gaussian smoothing. International Journal of Computational Intelligence in Bioinformatics and Systems Biology 1(3):316–331
  * Yao et al [2021] Yao L, Zhang J, Liu J, et al (2021) A deep learning-based system for bile duct annotation and station recognition in linear endoscopic ultrasound. EBioMedicine 65:103,238
  * Zhang et al [2020] Zhang J, Zhu L, Yao L, et al (2020) Deep learning–based pancreas segmentation and station recognition system in eus: Development and validation of a useful training tool (with video). Gastrointestinal endoscopy 92(4):874–885


$12700824$     $0.4845$     $49.61259375$     $19.485$     $1.75075$     $49.6338295$     $1.516125$     $0.0045125$     $0.0799125$

9|c|Average number of free channels ($F$) is $39$
$12448128$     $1.455$     $48.6255$         $137.89$      $27.355$       $48.790745$       $19.87375$     $0.072$        $1.10475$

$12459888$     $1.429$     $48.6714375$      $131.225$     $27.9175$      $48.83058$        $21.2225$      $0.07325$      $1.158875$

$12486256$     $1.428$     $48.7744375$      $135.71$      $27.91375$     $48.93806125$     $21.02$        $0.073125$     $1.1525$

$12596912$     $1.463$     $49.2066875$      $138.17$      $28.785$       $49.3736425$      $20.625$       $0.07475$      $1.146125$

$12650768$     $1.465$     $49.4170625$      $133.615$     $27.2375$      $49.577915$       $20.34$        $0.070375$     $1.112375$

$12710512$     $1.4$       $49.6504375$      $138.69$      $26.57125$     $49.81569875$     $20.06125$     $0.0685$       $1.091125$

$12828824$     $1.396$     $50.11259375$     $136.1$       $25.3825$      $50.27407625$     $18.91625$     $0.06475$      $1.02425$

$12903144$     $1.419$     $50.40290625$     $134.345$     $26.0275$      $50.56327875$     $18.91625$     $0.066$        $1.033$

$12927384$     $1.501$     $50.49759375$     $136.59$      $27.545$       $50.66172875$     $19.88125$     $0.06975$      $1.085375$

$12996424$     $1.425$     $50.76728125$     $136.87$      $27.82375$     $50.931975$       $20.9325$      $0.070125$     $1.123125$

$12700824$     $1.4381$     $49.61259375$     $135.9205$     $27.255875$     $49.775770125$     $20.178875$     $0.0702625$     $1.10315$

Average Delay and Delay Jitter for different real-life Text files

1.75cmData File Size (in bits)     1.75cmAverage Number of Channel Deallocation by PUs     1.75cmIdeal Transmission Time (IT) (in sec)     1.75cmInitial Channel Allocation Time (ICA) (in msec)     1.75cmChannel Reallocation Time (CR) (in msec)     2cmActual Transmission Time (AT) (in sec) (AT=IT+ICA+CR)     1.75cmMaximum Jitter (in msec)     1.75cmMean Jitter (in msec)     1.75cmStandard Deviation In Jitter (in msec)

9|c|Average number of free channels ($F$) is $697$
$25848$       $0$         $0.2019375$      $3.9925$     $0$           $0.20593$         $0$           $0$                $0$

$29008$       $0$         $0.226625$       $4.1225$     $0$           $0.2307475$       $0$           $0$                $0$

$30752$       $0$         $0.24025$        $4.035$      $0$           $0.244285$        $0$           $0$                $0$

$112920$      $0.004$     $0.8821875$      $4.11$       $0.00625$     $0.88630375$      $0.00625$     $0.0008928625$     $0.002375$

$322448$      $0.009$     $2.519125$       $4.0375$     $0.015$       $2.5231775$       $0.015$       $0.00075$          $0.003375$

$453664$      $0.014$     $3.54425$        $3.915$      $0.0275$      $3.5481925$       $0.0275$      $0.0009821375$     $0.00525$

$564880$      $0.013$     $4.413125$       $3.9925$     $0.02625$     $4.41714375$      $0.02625$     $0.00075$          $0.004375$

$703728$      $0.018$     $5.497875$       $4.0775$     $0.0275$      $5.50198$         $0.0275$      $0.0006395375$     $0.00425$

$1284416$     $0.029$     $10.0345$        $3.9675$     $0.05125$     $10.03851875$     $0.05125$     $0.0006487375$     $0.00575$

$3098904$     $0.06$      $24.2101875$     $3.9575$     $0.10625$     $24.21425125$     $0.105$       $0.0005592125$     $0.007625$

$662656.8$     $0.0147$     $5.17700625$     $4.02075$     $0.026$     $5.181053$     $0.025875$     $0.0005222488$     $0.0033$

9|c|Average number of free channels ($F$) is $510$
$25848$       $0.003$     $0.2019375$      $5.415$      $0.0075$      $0.20736$        $0.0075$      $0.00375$          $0.00525$

$29008$       $0$         $0.226625$       $5.565$      $0$           $0.23219$        $0$           $0$                $0$

$30752$       $0$         $0.24025$        $5.495$      $0$           $0.245745$       $0$           $0$                $0$

$112920$      $0.003$     $0.8821875$      $5.405$      $0.00625$     $0.88759875$     $0.00625$     $0.0008928625$     $0.002375$

$322448$      $0.008$     $2.519125$       $5.505$      $0.02$        $2.52465$        $0.02$        $0.001$            $0.0045$

$453664$      $0.016$     $3.54425$        $5.4275$     $0.03625$     $3.54971375$     $0.03625$     $0.00125$          $0.006875$

$564880$      $0.012$     $4.413125$       $5.6$        $0.02375$     $4.41874875$     $0.02375$     $0.000678575$      $0.004$

$703728$      $0.03$      $5.497875$       $5.4475$     $0.06375$     $5.50338625$     $0.06$        $0.0015$           $0.009375$

$1284416$     $0.04$      $10.0345$        $5.5$        $0.0975$      $10.0400975$     $0.09375$     $0.001234175$      $0.01075$

$3098904$     $0.076$     $24.2101875$     $5.3725$     $0.1825$      $24.2157425$     $0.175$       $0.000960525$      $0.012875$

$662656.8$     $0.0188$     $5.17700625$     $5.47325$     $0.04375$     $5.18252325$     $0.04225$     $0.0011266138$     $0.0056$

9|c|Average number of free channels ($F$) is $285$
$25848$       $0.002$     $0.2019375$      $9.1725$     $0.00375$     $0.21111375$      $0.00375$     $0.001875$     $0.002625$

$29008$       $0.001$     $0.226625$       $9.54$       $0.00375$     $0.23616875$      $0.00375$     $0.001875$     $0.002625$

$30752$       $0.001$     $0.24025$        $9.6875$     $0.0025$      $0.24994$         $0.0025$      $0.00125$      $0.00175$

$112920$      $0.006$     $0.8821875$      $9.2875$     $0.02125$     $0.89149625$      $0.02125$     $0.003$        $0.008$

$322448$      $0.018$     $2.519125$       $9.3875$     $0.0875$      $2.5286$          $0.0875$      $0.004375$     $0.019625$

$453664$      $0.028$     $3.54425$        $9.2975$     $0.08875$     $3.55363625$      $0.0875$      $0.003125$     $0.0165$

$564880$      $0.036$     $4.413125$       $9.7025$     $0.17125$     $4.42299875$      $0.17125$     $0.004875$     $0.029$

$703728$      $0.038$     $5.497875$       $9.495$      $0.11125$     $5.50748125$      $0.11125$     $0.002625$     $0.017$

$1284416$     $0.057$     $10.0345$        $9.5725$     $0.28125$     $10.04435375$     $0.2775$      $0.003625$     $0.0315$

$3098904$     $0.135$     $24.2101875$     $9.2675$     $0.54625$     $24.22000125$     $0.5175$      $0.002875$     $0.038125$

$662656.8$     $0.0322$     $5.17700625$     $9.441$     $0.13175$     $5.186579$     $0.128375$     $0.00295$     $0.016675$

9|c|Average number of free channels ($F$) is $39$
$25848$       $0.008$     $0.2019375$      $65.48$       $0.22875$     $0.26764625$     $0.22875$     $0.114375$     $0.16175$

$29008$       $0.011$     $0.226625$       $66.48$       $0.3425$      $0.2934475$      $0.3425$      $0.17125$      $0.242125$

$30752$       $0.01$      $0.24025$        $63.5325$     $0.37625$     $0.30415875$     $0.37625$     $0.188125$     $0.266$

$112920$      $0.044$     $0.8821875$      $64.85$       $1.05375$     $0.94809125$     $1.05375$     $0.1505$       $0.39825$

$322448$      $0.13$      $2.519125$       $66.66$       $4.40625$     $2.59019125$     $4.33375$     $0.220375$     $0.973$

$453664$      $0.162$     $3.54425$        $65.68$       $5.44375$     $3.61537375$     $5.19125$     $0.194375$     $0.994125$

$564880$      $0.204$     $4.413125$       $66.3075$     $6.0925$      $4.485525$       $5.8725$      $0.174125$     $1.001125$

$703728$      $0.215$     $5.497875$       $66.3075$     $6.04$        $5.5702225$      $5.74125$     $0.1405$       $0.884625$

$1284416$     $0.358$     $10.0345$        $66.82$       $9.1375$      $10.1104575$     $8.3625$      $0.115625$     $0.962125$

$3098904$     $0.508$     $24.2101875$     $66.67$       $10.94$       $24.2877975$     $9.53625$     $0.057625$     $0.716625$

$662656.8$     $0.165$     $5.17700625$     $65.87875$     $4.406125$     $5.247291125$     $4.103875$     $0.1526875$     $0.659975$

[1]
P. B. Galvin A. Silberschatz and G. Gagne.
Operating System Concepts.
John Wiley & Sons, New Delhi, eighth edition, 2009.

[2]
D. P. Agrawal and Qing-An Zeng.
Chapter 8 : Traffic channel allocation.
Online Lecture Notes, 2003.

[3]
E. Ahmed, A. Gani, S. Abolfazli, L.J. Yao, and S.U. Khan.
Channel assignment algorithms in cognitive radio networks: Taxonomy,
open issues, and challenges.
Communications Surveys Tutorials, IEEE, PP(99):1–1, 2014.

[4]
Ian F Akyildiz, Won-Yeol Lee, and Kaushik R Chowdhury.
Crahns: Cognitive radio ad hoc networks.
AD hoc networks, 7(5):810–836, 2009.

[5]
I.F. Akyildiz, Won-Yeol Lee, M.C. Vuran, and S. Mohanty.
A survey on spectrum management in cognitive radio networks.
Communications Magazine, IEEE, 46(4):40 –48, april 2008.

[6]
I.F. Akyildiz, W.Y. Lee, M.C. Vuran, and S. Mohanty.
Next generation/dynamic spectrum access/cognitive radio wireless
networks: a survey.
Computer Networks, 50(13):2127–2159, 2006.

[7]
Suzan Bayhan and Fatih Alagöz.
Distributed channel selection in crahns: A non-selfish scheme for
mitigating spectrum fragmentation.
Ad Hoc Networks, 10(5):774 – 788, 2012.
Special Issue on Cognitive Radio Ad Hoc Networks.

[8]
Suzan Bayhan and Fatih Alagöz.
A markovian approach for best-fit channel selection in cognitive
radio networks.
Ad Hoc Networks, 12(0):165 – 177, 2014.
(1) {SCEnarios} for ad hoc Network Evaluation Studies(SCENES); (2)
Modeling, Analysis and Simulation of Wireless and Mobile Systems & (3)
Algorithmic Aspects of Ad-Hoc Networks.

[9]
S. Bokharaiee, H.H. Nguyen, and E. Shwedyk.
Blind spectrum sensing for ofdm-based cognitive radio systems.
Vehicular Technology, IEEE Transactions on, 60(3):858–871,

[10]
M.M. Buddhikot.
Understanding dynamic spectrum access: Models,taxonomy and
In New Frontiers in Dynamic Spectrum Access Networks, 2007.
DySPAN 2007. 2nd IEEE International Symposium on, pages 649 –663, april

[11]
University of colorado boulder, Nov 2015.

[12]
ET FCC.
Docket no. 02-135 spectrum policy task force report, November 2002.

[13]
ET FCC.
Docket no 03-222 notice of proposed rule making and order, December

[14]
ET FCC.
Docket no 03-237 notice of inquiry and notice of proposed rulemaking,
November 2003.

[15]
A Behrouz Forouzan.
Data Communications & Networking (sie).
Tata McGraw-Hill Education, 2006.

[16]
Google images, Nov 2015.

[17]
Mp3 converted from youtube video, Nov 2015.

[18]
Youtube video, Nov 2015.

[19]
Zhao Hang-sheng and Wang Fan.
Qos provisioning single-channel opportunistic spectrum access
strategy in cognitive radio networks.
In High Performance Computing and Communications (HPCC), 2011
IEEE 13th International Conference on, pages 966–970, Sept 2011.

[20]
F.K. Jondral.
Software-defined radio: basics and evolution to cognitive radio.
EURASIP J. on Wireless Comm. and Networking, 2005(3):275–283,

[21]
A.S. Kamil and I. Khider.
Open research issues in cognitive radio.
architecture, 12:13, 2008.

[22]
Y.R. Kondareddy, P. Agrawal, and K. Sivalingam.
Cognitive radio network setup without a common control channel.
In IEEE MILCOM 2008, pages 1 –6, nov. 2008.

[23]
S. Krishnamurthy, M. Thoppian, S. Venkatesan, and R. Prakash.
Control channel based mac-layer configuration, routing and situation
awareness for cognitive radio networks.
In IEEE MILCOM 2005, pages 455–460, 2005.

[24]
H. Kushwaha, Yiping Xing, R. Chandramouli, and H. Heffes.
Reliable multimedia transmission over cognitive radio networks using
fountain codes.
Proceedings of the IEEE, 96(1):155–165, Jan 2008.

[25]
H. Mahmoud, T. Yucek, and H. Arslan.
Ofdm for cognitive radio: merits and challenges.
Wireless Communications, IEEE, 16(2):6–15, 2009.

[26]
A. Masri, C.F. Chiasserini, C. Casetti, and A. Perotti.
Common control channel allocation in cognitive radio networks through
uwb multihop communications.
In 1st Nordic Workshop. on Cross-Layer Optimization in Wireless
Networks, 2010.

[27]
MP Mishra and PC Saxena.
Survey of channel allocation algorithms research for cellular
Int J Netw Commun, 2(5):75–104, 2012.

[28]
III Mitola, J.
Cognitive radio for flexible mobile multimedia communications.
In Mobile Multimedia Communications, 1999. (MoMuC '99) 1999 IEEE
International Workshop on, pages 3 –10, 1999.

[29]
Hsien-Po Shiang and Mihaela van der Schaar.
Dynamic channel selection for multi-user video streaming over
cognitive radio networks.
In 15th IEEE International Conference on Image Processing,,
pages 2316–2319, Oct 2008.

[30]
Andrew S Tanenbaum.
Computer networks prentice hall.
New Jersey, 2006.

[31]
Shaowei Wang.
Efficient resource allocation algorithm for cognitive ofdm systems.
Communications Letters, IEEE, 14(8):725–727, 2010.

[32]
T.A. Weiss and F.K. Jondral.
Spectrum pooling: an innovative strategy for the enhancement of
spectrum efficiency.
Communications Magazine, IEEE, 42(3):S8–14, Mar 2004.

[33]
Chunsheng Xin and Xiaojun Cao.
A cognitive radio network architecture without control channel.
In Proc. IEEE GLOBECOM 2009, pages 1–6, Nov 2009.

[34]
Lei Yang, Wei Hou, Lili Cao, Ben Y Zhao, and Haitao Zheng.
Supporting demanding wireless applications with frequency-agile
In NSDI, pages 65–80, 2010.

[35]
Yuan Yuan, Paramvir Bahl, Ranveer Chandra, Thomas Moscibroda, and Yunnan Wu.
Allocating dynamic time-spectrum blocks in cognitive radio networks.
In Proc. MobiHoc, pages 130–139, 2007.

[36]
Shuai Zhang, Yanhui Lu, Fang Wang, and Xiaomin Mu.
Spectrum access strategies based on the ratio of throughput to
spectrum consumption in the context of spectrum fragmentation in distributed
cognitive radio networks.
In Communication Systems (ICCS), 2014 IEEE International
Conference on, pages 442–446, Nov 2014.

[37]
Dayong Zhou and M. Zukerman.
Performance and efficiency evaluation of channel allocation schemes
for hscsd in gsm.
In Proc. IEEE GLOBECOM, volume 1B, pages 1084–1088 vol. 1b,

[38]
M. Zukerman.
Bandwidth allocation for bursty isochronous traffic in a hybrid
switching system.
Communications, IEEE Transactions on, 37(12):1367–1371, Dec