# Open-Source, Cost-Aware Kinematically Feasible Planning
for Mobile and Surface Robotics

Steve Macenski  Open Navigation LLC
<EMAIL_ADDRESS>Matthew Booker  Carnegie Mellon University
<EMAIL_ADDRESS>Joshua Wallace  Locus Robotics
<EMAIL_ADDRESS>
###### Abstract

This paper introduces the Smac Planner, an openly available search-based
planning framework with multiple algorithm implementations including 2D-A*,
Hybrid-A*, and State Lattice planners. This work is motivated by the lack of
performant and available feasible planners for mobile and surface robotics
research.

This paper contains three main contributions. First, it briefly describes a
minimal open-source software framework where search-based planners may be
easily added. Further, this paper characterizes new variations on the feasible
planners - dubbed Cost-Aware \- specific to mobile roboticist’s needs. This
fills the gap of missing kinematically feasible implementations suitable for
academic, extension, and deployed use. Finally, we provide baseline
benchmarking against other standard planning frameworks.

Smac Planner has further significance by becoming the standard open-source
planning system within ROS 2’s Nav2 framework which powers thousands of robots
in research and industry.

##  I Introduction

Path planning for mobile platforms has undertaken significant development in
recent years. The classical infeasible or holonomic techniques have made way
for a more mature class of algorithms. These can provide additional features
such as guaranteeing drivability, respecting kinodynamic constraints, and
operating in advanced environmental representations.

The DARPA Grand and Urban Challenges created innovation in planning to support
Ackermann models [1]. These competitions can be accredited with the creation
or refinement of techniques such as Hybrid-A* and State Lattice.

Due to Moore’s law, many of these techniques are now applicable to mobile and
surface marine robot systems which contain far more limited compute resources.
Search-based methods continue to be efficient and effective in low-dimensional
state spaces. However, mobile robots continue to rely on these infeasible
techniques such as Dikjstra’s Algorithm, Navigation Function, and numerous
forms of grid-based heuristic search (e.g. A*, D*) [2, 3]. For circular
differential-drive or holonomic robots, this class of algorithm is still
suitable for many uses.

Figure 1: Large, non-circular industrial robot using Nav2 uniquely enabled by
this open-source Framework.

However, the sector has evolved away from simple circular robots and towards
new industries outside of warehouses and research laboratories. There is a
growing use of Ackermann, large non-circular, and quadruped robots in new
applications such as last-mile delivery, automated platform lifts, and
construction monitoring. These robots cannot utilize single-point holonomic
planning techniques previously deployed and require, at minimum,
kinematically-feasible paths. Yet, common feasible algorithms have no quality
implementations accessible for academic use, are not integrated with any major
research framework (MRPT, ROS, etc), nor have variants unique to mobile
roboticist’s needs [4, 5].

Our work was motivated by this gap: to serve these new styles of robots and
its researchers. To fill this, we created the Smac Planner: a templated,
C++-based search-based planning framework designed to make creating search-
based planning algorithms simple. Within it, we propose three planners; 2D-A*,
Hybrid-A*, and State Lattice.

This paper contains three areas of contribution. It describes the high-level
design of the Smac Planner framework, enabling performant planners to be added
in approximately 200 lines of C++. More substantially, this work characterizes
several variations of the kinematically feasible planners (Hybrid-A*, State
Lattice) to build upon them for the requirements and conventions in mobile
robotics. Finally, we provide performance benchmarking and analysis of the
planners.

We provide high-quality implementations in Nav2 111https://github.com/ros-
planning/navigation2 for rapid benchmarking and evaluation [6]. This is in use
by over a dozen organizations worldwide in research, warehouses, outdoors, and
in surface marine applications; including the industrial warehouse robot in
Fig 1. Its presence has enabled Nav2 to be integrated with a broader class of
robots including Ackermann, legged, and large non-circular robots,
accelerating research on these platforms [7].

##  II Related Work

Systems that are not able to be modeled as circular or holonomic require more
developed planning techniques for robust operation. Kinematically feasible and
kinodynamic planners model kinematics and/or dynamics to provide plans for
limited maneuverability drive-trains and arbitrarily shaped robots [8].
Popular search-based feasible planners include Hybrid-A* and State Lattice.
Hybrid-A* makes use of Ackermann curvature-constrained models such as Dubins
and Reeds-Shepp to search grid maps. Meanwhile State Lattice applies a set of
pre-generated motion primitives, known as the control set, to find its
neighbors for arbitrary motion models [9, 10]. As such, it is required to
generate the control set, usually offline. One such algorithm is described in
[10] and in later detail (Sec. V-C).

Sampling-based methods are common for kinodynamic planning, however they have
slow run-time performance over large distances in the presence of non-trivial
obstacles (e.g. long time horizons) [12]. Favorably, hybrid planning used in
navigation systems creates a distinction between global path planning and
local trajectory planning such that kinodynamic planning is leveraged only in
the local time horizon. Generally, global planning need only apply kinematic
constraints [13].

Many planning frameworks have been proposed, notable among them with staying
power are the Open Motion Planning Library (OMPL) and the Search-Based
Planning Library (SBPL) [14, 15]. OMPL specializes in sampling-based
techniques with implementations of common methods, including PRM, RRT and its
variants. While OMPL can sample using mobile robot models, in experiments, it
provided plans with irregular run-times often an order of magnitude (or more)
slower than search-based methods for the low-dimensional and non-trivially
occupied state spaces. While technically drivable, it suffers from significant
unnecessary sudden turns - including long segments of reversing near obstacles
even when both are explicitly penalized.

SBPL instead focuses on search-based planning with several heuristic search
options such as Anytime Repairing A* and D* Lite [17]. It uses motion
primitives to search, but the primitive generators within SBPL are manually-
engineered and do not create principled minimum control sets, limiting
maneuverability. SBPL also only provides naive distance-based heuristics,
resulting in acutely slow planning times in structured, realistic
environments. The current state of this framework is more suitable for local
trajectory planning - as recent publications would also suggest [15].

Cost maps or grids are a common environmental model which balances fidelity
with efficiency to aggregate information about a scene from various sources
[16]. While they evolved from early probability occupancy grids, the modern
cost map uses values to represent relative costs-to-go rather than likelihoods
of occupancy. Though, there are special values reserved for occupancy to
provide hard collision constraints. These cost grids are common to mobile
robotic navigation systems as the model of the world for which planning and
control takes place.

Prior to this work, the primary planner in the Nav2 framework is the NavFn
planner [6]. It implements a Navigation Function planner that uses cost grid
information in its search and collision checks based on point-costs. Is is
unable to navigate a highly asymmetric robot through narrow spaces or provide
any feasibility guarantees. While NavFn is efficient, it lacks the necessary
features for Ackermann, legged, or non-circular robot systems - representing a
significant gap. All path and trajectory planners in the framework utilize the
cost grid information to generate plans.

##  III Smac Planner Framework

The Smac Planner consists of three fundamental layers - a templated A* for
search, node planner templates, and an integration layer. A* was selected due
to its ease of integration with arbitrary graph types, search neighborhoods,
and heuristic computations needed for a variety of planners. The A* is
templated by the node planner type, NodeT in Fig. 2, which implements the
planner-specifics. Different templates use the decoupled and abstracted A*
algorithm to create new planners with different characteristics, which
minimizes boilerplate. This abstraction decouples the cost functions and
neighbor selection policies of a particular planner from the shared A*
algorithm. While this framework is not particularly novel, it is beneficial to
expand briefly on its design as a reference to the openly available
implementation also provided in 1.

Figure 2: Outline of the A* and key node template methods.

Node planner types serve two purposes - they contain the node’s state in the
graph and implement the planner. As a node, it stores the typical data such as
the collision state, if it has been queued or visited, and the accumulated
cost.

It also contains the planner-specific logic used to compute traversal costs,
heuristics, and expansion characteristics, shown in Fig. 2. The non-state data
associated with the planner are shared across all node instances using static
members to gain the performance of a standalone implementation. New planners
are created as NodeT’s to implement the appropriate search methodology.

The link to a robotics framework is created only at the highest level
integration layer using an instantiation of the A* algorithm with the
appropriate template. Thereby, this work may be applied to multiple robotics
frameworks easily and has been done so.

##  IV Variations for Mobile Robots

Mobile robots generally have a different operating context than autonomous
vehicles (AV) on roads. They have less structured interactions with other
agents in the environment - often without lanes, intersections, or speed
conventions. Additionally, their free-space planners must run on more
restricted compute and provide tunable behaviors for a broad range of
applications.

This section explores a few alterations applied to the feasible AV planners
motivated by the application to mobile robotics. Chief amongst them is the
Cost-Aware Obstacle Heuristic, giving the planner variants their namesake.

###  IV-A Cost-Aware Obstacle Heuristic

The Cost-Aware Obstacle Heuristic shepherds the expensive feasible search away
from obstacles, towards optimal path solutions, while respecting soft grid map
constraints. It performs a 2D-grid search and caches the accumulated path
costs to use as the heuristic. It reduces planning time by avoiding dead-ends,
U-turns, and wrong directions using the obstacle occupancy information in the
grid map. Cost-aware refers to utilizing all cost map values rather than only
the binary obstacle information, allowing behavioral constraints placed in the
cost grid to be accounted for during planning.

The heuristic is computed from the goal pose so that it can be expanded
dynamically and cached. It uses a dynamic programming implementation of a
Differential A* so that preexisting queued nodes may be re-prioritized and
heuristic search restarted during the feasible planner’s execution [3].

The 2D search’s traversal cost uses the distance traveled and an additional
weighted cost-proportional term, in Eq. 1.

$cost_{traversal}=dist*(1+\frac{w_{cost}*cost_{i,j}}{cost_{max}})\\\ $ (1)

where $w_{cost}$ is the path cost penalty, $cost_{i,j}$ is the cost of moving
from cell i to cell j and $cost_{max}$ is the maximum value in the cost map.
Normalizing the grid cost scales it proportionally to the travel distance,
while the penalty weight tunes the trade off between path cost and path
length. This can be interpreted as a percentage increase in cost-to-travel
based on the grid cost. The planners’ traversal costs are formulated the same
as Eq. 1 using cost-weighted distances. Thereby this 2D cost-weighted search
heuristic is admissible.

This heuristic has an equally admissible analog in classical Hybrid-A* and
others which takes into account only lethal occupancy information [9, 18].
They instead use an expensive path smoothing algorithm to maximize the
distance from all obstacles irrespective of context using a Voronoi diagram of
the map. The State Lattice technique has no such analogous heuristic.

Figure 3: Cost-Aware Obstacle Heuristic (left) vs Obstacle Heuristic (right).

Our method has several advantages for mobile robotics. Often, the behavior of
the mobile robot is created through costs and potential fields in a grid map.
These costs create desirable behaviors, such as punishing traversal in front
of other dynamic agents, inflating obstacles to punish traversal close to
potential collisions, or rewarding exploration in certain regions over others.
In that way, robot systems define more complex behaviors and constraints in
the environment than simply maximizing distance from obstacles, which could be
effectively modelled in a Voronoi diagram.

Instead, it is favorable to use this contextual information in the search
itself for mobile robotics applications. Many behaviors are conventionally
conveyed in the grid map, not as least squares constraints for post-search
smoothing. Further, behaviors similar to those in Fig. 3 explore an entirely
separate solution space and would be challenging to obtain reliably with only
post-processing.

Fig. 3 shows an illustrative example where higher costs are applied to a
region to dissuade travel. This soft constraint may be introduced due to
increased risk or danger in that area (red). On the left, the Cost-Aware
Obstacle Heuristic uses this constraint and finds a longer path, going around
the danger zone to achieve a lower path cost. On the right, the binary
Obstacle Heuristic ignores this constraint and marches through the region for
a slight improvement in path length. Similarly, a potential field penalty was
added to the walls. Our heuristic uniquely steered towards the goal and
maintained acceptably navigable distances from obstacles.

As this heuristic has cost-awareness, it curtails the strict requirement for
expensive path smoothing. Indeed, using an aggressive smoothing algorithm
could easily erase the behaviors and constraints considered during search.

###  IV-B Traversal Penalty Functions

In addition to constraints in the grid map from external factors, it is
beneficial to constrain the search for path quality. Thus, we define a set of
trivial penalty functions in the traversal function to influence the planners’
search behavior. This, in conjunction with the Cost-Awareness, removes the
requirement for expensive path post-processing for smoothness when well tuned.

The first of these is the cost penalty discussed in Sec. IV-A. This penalty is
unique in that it is applied in both the cost-aware heuristic and planners’
traversal functions for admissibility. Thus, we can study the impact of
varying heuristic behavior by analyzing the Smac 2D-A* planner in isolation
with varying cost penalties, shown in Fig 5. Eq. 1 indicates larger values of
this penalty will punish expanding search into higher cost areas at the
expense of longer paths. Set to $0.0$, the cost-awareness is disabled and the
optimal path skirts the walls. With modest weights applied ($1.0$-$4.0$), the
path displays a more rational approach into the narrow corridor on the lower
right and favors open spaces to varying degrees. The higher the value in the
$1.0$-$4.0$ range, the longer the path in attempt to minimize proximity to the
walls. Sufficiently high values ($5.0$+) punishes narrow spaces severely - as
$8.0$ explores another route solution entirely.

The next two penalty functions are the non-straight and change penalties.
These two penalties correct for any wobbling or unnecessary turning behavior.
The non-straight penalty applies a cost to any not straight motion primitive.
The change penalty is complementary - it punishes changing turning directions
($\Delta\alpha$) to commit to an action. A small penalty helps create smooth
and reliable behavior ($2$-$5\%$ yields good results) but larger values can
create straight-line and rigid behavior which may be beneficial in some
settings. These penalties are multiplied with the $cost_{trav.}$ in Eq. 1 to
find the final travel cost. When both changing direction and non-straight, the
penalties are summed, Eq. 2.

$cost_{travel}=$ $\begin{cases}cost_{trav.}&\alpha_{i}=0\\\
(1+w_{nonstrt})\hskip
2.84526ptcost_{trav.}&sgn(\alpha_{i})=sgn(\alpha_{i-1})\\\
(1+w_{nonstrt}+w_{ch})\hskip 2.84526ptcost_{trav.}&sgn(\alpha_{i})\neq
sgn(\alpha_{i-1})\\\ \end{cases}$ (2)

where $\alpha_{i}$ is the change in turning direction, $w_{nonstrt}$ is the
non-straight path penalty, and $w_{ch}$ is the penalty for flipping directions
from turning left to right or vice-versa.

The final traversal penalty is reverse. This penalizes any motion primitive in
the reverse direction to dissuade backing up. A value of $2.1$ allows
reversing to turn around when given a goal in the opposite heading direction,
but solemnly used otherwise. This is applied after all other penalties.

###  IV-C Search Resolution

The grid map’s primary purpose is to provide a search space for planning. An
architect must make a determination about the resolution of a grid map based
on the fidelity of sensor data, type of environment (such as confined spaces),
and compute available.

Originally, Hybrid-A* was introduced using a 0.5m search grid on a 0.15m
obstacle map to achieve 50-300ms performance - effectively reducing the chosen
search resolution and increasing the primitive lengths by $\sim$3.3x. Then, it
smooths and upsamples to create a plan at the original resolution. Fig 4
analyzes paths generated using Hybrid search on a Dubin’s model, with
primitives proportional to 1x and 3x the grid resolution with otherwise
identical settings. The smaller primitives at grid resolution create more
stable search behavior while the 3x effectively requires expensive smoothing
to recover an acceptable path. Further, lengthy primitives can make planning
in narrow passages difficult.

Figure 4: Hybrid-A* primitive size vs unsmoothed path quality.

Thus, all of the Smac Planners use the grid resolution to adjust their motion
primitives. Naturally, the 2D-A* uses its neighboring cells exactly. Smac’s
Hybrid-A* planner dynamically generates a set of primitives dependent on the
grid resolution to guarantee that its length will leave the current grid cell
and is at least one angular quantization away (for turning motions). The State
Lattice implementation computes a minimum control set offline. Each of its
primitives are non-trivially longer than the grid resolution, so they are
distretized by the grid map resolution for collision detection and path
reconstruction.

##  V Path Planners

The planners are created as node planner types (NodeT) on A*. The framework
provides three implementations; 2D-A*, Hybrid-A*, and State Lattice. Each are
implemented in under 300 lines of planner-specific code with an average of
188; additional templates can be easily created.

###  V-A Cost-Aware 2D-A* Planner

A Cost-Aware 2D-A* algorithm is included to provide baseline functionality for
circular robots. It is similar to the heuristic search from Sec. IV-A that
provides the feasible planners with its Cost-Aware search guide. In fact, it
has identical traversal and heuristic properties to the Cost-Aware Obstacle
Heuristic, but is structured as an end-user 2D planner that does not cache its
heuristic field to save compute resources. Hence it shares the Cost-Aware
moniker.

Figure 5: Cost-Aware 2D-A* with variable cost penalties.

It uses Moore (8) connected spaces for search with an L2 distance heuristic,
which is admissible for any possible value of cost. Fig. 5 shows the Cost-
Aware 2D-A* in a warehouse environment finding paths approximately 85m in
length.

###  V-B Cost-Aware Hybrid-A* Planner

The Cost-Aware Hybrid-A* planner provides feasible paths using a Dubins or
Reeds-Shepp motion model constrained by a minimum turning radius for ackermann
vehicles. It stores the continuous coordinates of the expanded nodes to create
a drivable path. This planner typically runs in 20-300ms across a warehouse
(approx. 6500$m^{2}$), dependent on the space’s complexity. Speeds as fast as
5ms are displayed when the heuristic is cached for replanning (optional).

It uses two admissible heuristics, of which it takes the maximum. The first is
the Cost-Aware Obstacle heuristic previously discussed in Sec. IV-A which
steers the expansion down useful directions with particular behaviors. The
second is the same as in the original paper’s ‘non-holonomic without
obstacles’ admissible distance heuristic, using a pre-computed lookup table of
Dubin or Reeds-Shepp distances from each grid cell and angular quantization
around the goal pose assuming no obstacles [9]. The lookup table uses a
technique introduced by State Lattice- storing half the window and using the
model’s symmetry to mirror poses [10]. This halves the memory usage for the
same sized heuristic lookup window and is newly applied to the Hybrid-A*
technique. This heuristic prunes approach headings to the goal that would be
difficult to achieve from the drivetrain’s constraints.

The planner’s traversal function is that from Sec. IV-B. Also as previously
discussed, the motion primitives are dynamically computed proportional to the
grid map resolution. The minimum required angular change ($\Delta\theta$) is
determined by the angle of a chord of at least $\sqrt{2}$ cells long from a
circle with minimum turning radius, $r_{min}$, Eq. 3. $\Delta\theta$ is
converted into angle quantization coordinates ($\Delta\theta_{q}$); the
position changes are determined geometrically. Each angle has its primitives
pre-computed and translated to avoid repeated trigonometric operations.

$\begin{gathered}\Delta\theta\geq 2\sin^{-1}({\dfrac{\sqrt{2}}{2r_{min}}})\\\
\Delta x=r_{min}\sin\Delta\theta_{q}\\\ \Delta
y=r_{min}(1-\cos\Delta\theta_{q})\\\ \end{gathered}$ (3)

Analytic expansion is employed to determine an obstacle-agnostic minimal path
to the goal from a currently expanded node. It is attempted at a higher
frequency in closer proximity to the goal using the heuristic as a measure.
New to this variant, the analytic path is only used if it is sufficiently
close to the goal. This continues the search process while still far from the
goal to account for the map constraints. While using the Reeds-Shepp model,
the maximum length also prevents returning auto-generated, long, backward
paths unless on approach to the goal, where it may be strictly required. The
expansion is finally utilized if collision-free.

Collision checking is performed on each end pose. Like the cached primitives,
the footprint is pre-rotated, then translated at run-time. When a potential
field is applied, an optimization may be used. We find the cost of the center
of the robot at which the footprint is possibly in collision. Then, a single
point check can indicate the pose’s feasibility in some cases. Because the
heuristic steers away from obstacles, the proportion of full footprint checks
is small.

###  V-C Cost-Aware State Lattice Planner

The Cost-Aware State Lattice planner searches using a set of motion primitives
generated using an arbitrary motion model’s constraints (e.g. ackermann,
differential, omnidirectional). The run-time of this planner is similar to
Cost-Aware Hybrid-A* at 10-350ms in the same environment. With the heuristics
cached for replanning, it approaches 1-3ms. Motion primitives may be created
using the provided generator, which builds principled minimum control sets
discussed below. However, other techniques may be substituted.

Similar to the previous section, this planner uses the same admissible
heuristics, continuous cell coordinates, pre-computation optimizations, and
traversal functions as the Cost-Aware Hybrid-A*. The analytic expansion and
cost-aware heuristic are both new to the state lattice planner regime and
significantly decrease planning times [10]. Differently, this planner’s
primitives are discretized by the grid resolution to construct dense paths and
evaluate intermediary poses for collision.

####  V-C1 Generating the Minimum Control Set

The minimal control set is the smallest set of motions needed to represent the
state lattice of a motion model. The method utilized removes pairs of states
with similar motions which can be considered redundant, condensing the state
lattice to its minimal set [11, 10]. A minimal control set generated using our
system is given in Fig. 6. In this figure, each of 16 angular quantizations
have between 3-5 unique motion primitives.

Figure 6: A control set using a 5cm grid resolution and a 1m turning radius
with 16 angle quantizations.

An important aspect in creating the minimal control set is path decomposition.
A path is decomposable if it can be replaced by the concatenation of paths in
the control set. This is evaluated by checking if the composed constituent
path has an end configuration within half a grid cell or angle orientation
quantization of the original. Thus, the minimal control set is what remains
when all decomposable paths are removed; yielding the minimum spanning
motions.

Generating the minimal control set is achieved by iterating over all end poses
and start headings from the origin to find feasible trajectories. If one
exists, it is checked for decomposition by trajectories in the current set and
added if unique. End points are selected based on a wavefront that moves away
from the origin. This process generates shorter paths first, leading to a
higher likelihood of subsequent paths being decomposable. Eventually, all
candidate trajectories in a given wavefront will be decomposable and add no
new information to the control set. The number of fully-represented
consecutive wavefronts is the termination condition.

One consideration when generating the set lies in the discretization of
headings. A uniform discretization leads to the generation of minimal control
sets with excessively long, non-straight paths [11]. Instead, center positions
of cells in the wavefront are used to derive non-uniform heading angles. In
doing so, useful motions are significantly shorter and decompose to fewer
consequential trajectories- creating a more tractable control set.

####  V-C2 Generating Trajectories in the Set

To generate control sets, we must generate candidate feasible trajectories to
any given state. We propose an intuitive method for trajectory generation that
guarantees minimal curvature, Fig. 7. Given a start configuration
$(x_{i},y_{i},\theta_{i})$ and an end configuration
$(x_{f},y_{f},\theta_{f})$, we define two lines, $l_{1}$ which passes through
the start configuration and $l_{2}$ which passes through the end
configuration. Then, we compute their intersection point, I, Eq. V-C2.

$\displaystyle a=$
$\displaystyle\frac{y_{f}-x_{f}\tan\theta_{f}-y_{i}+x_{i}\tan\theta_{i}}{\tan\theta_{i}-\tan\theta_{f}}$
$\displaystyle I=$ $\displaystyle[a,l_{1}(a)]$ (4)

Next, we calculate the end points of the trajectories arc by Eq. V-C2.
$p^{arc}$ will lie on $l_{1}$ and $q^{arc}$ will lie on $l_{2}$. The end
points are calculated using distance $d$, the minimum of either I to p or I to
q, where $p=(x_{i},y_{i})$ and $q=(x_{f},y_{f})$. By nature of $d$, either
$p^{arc}$ will be coincident with $p$, or $q^{arc}$ will be coincident with
$q$.

$\displaystyle d=$ $\displaystyle\min(\lVert I-p\rVert,\lVert I-q\rVert)$
$\displaystyle p^{arc}=$ $\displaystyle d\>\overrightarrow{p-I}+I$
$\displaystyle q^{arc}=$ $\displaystyle d\>\overrightarrow{q-I}+I$ (5)

Finally, we find the center point $C$ and radius $r$ of the circle whose arc
joins $p^{arc}$ at an angle of $\theta_{i}$, to $q^{arc}$ at $\theta_{f}$.
This is found via an intersection point of perpendicular bisectors of $l_{1}$
at $p^{arc}$ and $l_{2}$ at $q^{arc}$. This is the center since the distances
from $p^{arc}$ and $q^{arc}$ are $r$ (Eq. V-C2).

$\displaystyle l_{1}^{\perp}(x)=$ $\displaystyle
x\frac{-1}{\tan\theta_{i}}+p^{arc}_{y}+\frac{p^{arc}_{x}}{\tan\theta_{i}}$
$\displaystyle l_{2}^{\perp}(x)=$ $\displaystyle
x\frac{-1}{\tan\theta_{f}}+q^{arc}_{y}+\frac{q^{arc}_{x}}{\tan\theta_{f}}$
$\displaystyle b=$
$\displaystyle\frac{q^{arc}_{y}+\frac{q^{arc}_{x}}{\tan\theta_{f}}-p^{arc}_{y}-\frac{p^{arc}_{x}}{\tan\theta_{i}}}{\frac{-1}{\tan\theta_{i}}+\frac{1}{\tan\theta_{f}}}$
$\displaystyle C=$ $\displaystyle[b,l_{1}^{\perp}(b)]$ (6) $\displaystyle r=$
$\displaystyle\lVert C-p^{arc}\rVert$ (7)

A straight-line segment is added to the start or end configuration connected
to the arc to complete the trajectory. Using this method, the trajectory is
either an arc or contains an arc and a line segment, depending on the inputs,
in which guarantees that the trajectory contains the minimal curvature.

Figure 7: An example trajectory for start configuration (0,0,$\frac{\pi}{8}$)
to end configuration (2,2,$\frac{\pi}{2}$).

An additional step ensures that the minimal curvature is within the vehicle’s
physical constraint, $R$. The closer $q^{arc}$ and $p^{arc}$ are to I, the
smaller the radius of the path. This places a lower limit on $d$ for the
trajectory to be valid. Eq. V-C2 gives the expression for the minimum
distance.

$\displaystyle\psi=$ $\displaystyle\frac{\pi}{2}-|\theta_{f}-\theta_{i}|$
$\displaystyle d_{min}=$ $\displaystyle\frac{R}{\tan\frac{\psi}{2}}$ (8)

##  VI Experiments and Analysis

###  VI-A Random Environments

This experiment benchmarks the performance of each of the planners in the
framework along with baseline planners NavFn and SBPL using ARA*. Three random
maps of 10,000$m^{2}$ with varying occupancy of 10%, 15%, 20% were generated
[19]. The compute time and paths were recorded for each plan. All experiments
were run on an AMD Ryzen 5 5600X CPU with Ubuntu 20.04.

A set of 1,000 start and goal pairs were randomly generated. The pairings had
at least 3m of distance between them and were verified to be reachable. Fig. 8
shows all of the paths generated on the 20% occupancy map.

Each of the Smac Planners were run with complimentary parameters: $r_{min}$ of
0.4m, reversing was allowed, a cost penalty of $2.0$ and a $5\%$ non-straight
and change penalties. SBPL’s primitive set was created using its provided
generator with equivalent parameters: 16 heading discretizations, 5cm grid
resolution and 5 primitives per heading. The results are shown in Table I.

TABLE I: Random occupancy map results. Best results and those within 1% are
marked in bold.

Hybrid-A* State Lattice 2D-A* SBPL NavFn $10\%$ $t(ms)$ 39.07 $42.31$ $66.15$
$5,640$ $71.11$ $l_{path}(m)$ 51.41 51.40 50.96 $51.99$ $52.60$ $15\%$ $t(ms)$
40.73 $43.25$ $85.63$ $6,587$ $66.45$ $l_{path}(m)$ $51.10$ $51.15$ 50.45
$54.87$ $52.50$ $20\%$ $t(ms)$ 38.77 $39.40$ $88.82$ $6,633$ $61.02$
$l_{path}(m)$ $50.78$ $50.51$ 49.65 $53.68$ $52.25$

Figure 8: All paths generated for the 20% occupied map.

Across the experiments, the feasible Smac Planners had the most stable
behavior. Both Cost-Aware Hybrid-A* and State Lattice had run-times between
39-42ms, on average. In fact, they demonstrated less than 5% difference in
run-time performance and less than 0.2% difference in paths lengths. This
similar behavior only appears to differ to the extent that their available
primitive sets attempts to approximate the same motion. This can be attributed
to the use of the same traversal and heuristic functions for both of the
search models. Further, they were both markedly faster than Cost-Aware 2D-A*
(approx. 50%) and NavFn (approx. 38%).

Figure 9: Example paths in the 20% occupied map.

The Cost-Aware 2D-A* had reliably the shortest paths. However, the lengths of
the Cost-Aware Hybrid-A* and State Lattice paths were within 2.5% of it in the
experiments. Cost-Aware 2D-A* has slightly shorter paths due to its
8-connected search model preferring straight line motions, rather than
restricted turning maneuvers. A slight increase in length is to be expected to
achieve feasibility.

SBPL’s state lattice ARA* planner displayed two orders of magnitude slower
global planning than any other method tested - 150x slower than its Cost-Aware
State Lattice counterpart. This is due to the lack of modern heuristics and
optimizations. Much of the compute time is spent searching on approach to the
goal at the right heading. SBPL also has longer paths than any of the proposed
planners, as well as 2/3 of NavFn’s evaluations. SBPL paths regularly have
more sharp, unsmooth turning behaviors.

Fig. 9 shows a set of example paths from the experiment with randomized start
and goal poses. The Cost-Aware 2D-A* is rather mechanical with straight-line
segments. Smac’s Hybrid-A* and State Lattice planners have significant overlap
and feasibly navigated the maze retaining suitably navigable distances from
obstacles. The NavFn planner also created a relatively smooth motion
maintaining its distance from obstacles. However, it still lacks in guaranteed
feasibility and is slower.

Overall, the experiment shows that the feasible Smac Planners have run-times
well in excess to its 2D and SBPL counterparts. All of the Smac Planners
generate paths shorter than both of the baselines - and the kinematically-
feasible Smac Planners are notably faster. The behavior of Cost-Aware
Hybrid-A* and State Lattice in this experiment is largely equivalent. Thus,
the selection between them may depend more on the motion model most desired
rather than specific architectures.

###  VI-B Real-World Warehouse

In this experiment, we tested the three Smac Planners in a 363,000$ft^{2}$
warehouse using an Intel i7-8565U CPU with same parameters as in Sec. VI.
While this map was generated using a proprietary system, similar may be
obtained using open tools like SLAM Toolbox [20]. The three goal poses in Fig.
10 were selected to display planning characteristics in both major and minor
aisles distributed across the map.

Figure 10: A set of example paths of the Smac Planner in a warehouse serviced
by Locus Robotics (33,600$m^{2}$).

It is clear that the Cost-Aware 2D-A* plans are suboptimal for execution and
make unnatural movements when outside the main thoroughfares. The feasible
planners create high-quality paths even in challenging confined environments
in the minor aisles. They both display instinctive maneuvers to transition
from one aisle to another. Unsurprisingly, both Hybrid-A* and State Lattice in
Smac demonstrate similar behavior.

The feasible planners also show better run-time performance. The Smac Planner
computation times averaged over 10 trials showed that Hybrid-A* is the fastest
at 290ms, State Lattice at 473ms, and the 2D-A* at 1358ms in this
exceptionally large and dense warehouse environment.

##  VII Conclusion

We introduced the Smac Planner, a minimal framework for building performant
search-based planners. Within it, we described our variations to create a set
of kinematically-feasible path planners specific to the needs of modern mobile
roboticists. The Cost-Aware planners contain variations in their heuristics,
traversal functions, and high level architecture to maximize performance and
follow conventions common to the mobile robotics field. Experiments
benchmarked the new planners’ performance in sandbox and real-world, large-
scale situations.

This fills the gap of missing feasible planners within a robotics frameworks.
These planners are in use today by over a dozen organizations worldwide,
enabling car-like, legged, and large non-circular robots for the first time in
ROS. The software behind this work can be found here 1.

## Acknowledgements

We would like to thank Locus Robotics for their generous assistance in
providing the warehouse map. We would also like to thank RoboTech Vision for
their detailed assistance.

## References

  * [1] B. Paden, M. Čáp, S. Z. Yong, D. Yershov and E. Frazzoli. ”A Survey of Motion Planning and Control Techniques for Self-Driving Urban Vehicles,” in IEEE Transactions on Intelligent Vehicles, vol. 1, no. 1, pp. 33-55, March 2016.
  * [2] J.-C. Latombe. ”Robot motion planning”. Kluwer Academic Publishers, Dordrecht, Netherlands, 1991.
  * [3] K. Trovato. ”Differential A*: An adaptive search method illustrated with robot path planning for moving obstacles and goals, and an uncertain environment”. International Journal of Pattern Recognition and Artificial Intelligence 4.02, 1990.
  * [4] S. Macenski, T. Foote, B. Gerkey, C. Lalancette, W. Woodall. ”Robot Operating System 2: Design, architecture, and uses in the wild”. Science Robotics, 7(66), eabm6074, May 2022.
  * [5] J. Claraco. ”Development of scientific applications with the mobile robot programming toolkit”, 2008.
  * [6] S. Macenski and F. Martin and R. White and J. Ginés Clavero. ”The Marathon 2: A Navigation System”. 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS).
  * [7] S. Macenski, T. Moore, DV Lu, A. Merzlyakov, M. Ferguson, ”From the desks of ROS maintainers: A survey of modern & capable mobile robotics algorithms in the robot operating system 2”, Robotics and Autonomous Systems, 2023.
  * [8] E. Schmerling, and M. Pavone. ”Kinodynamic Planning”. Encyclopedia of Robotics, Springer, 2019.
  * [9] D. Dolgov, S. Thrun, M. Montemerlo and J. Diebel. ”Practical Search Techniques in Path Planning for Autonomous Driving”. AAAI Workshop - Technical Report, 2008.
  * [10] M. Pivtoraiko, R. Knepper, and A. Kelly. ”Optimal, smooth, nonholonomic mobile robot motion planning in state lattices”. Robotics Institute, Carnegie Mellon University, Pittsburgh, PA, Tech. Rep. CMU-RI-TR-07-15, 2007.
  * [11] M. Pivtoraiko and A. Kelly, ”Generating near minimal spanning control sets for constrained motion planning in discrete state spaces,” 2005 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2005, pp. 3231-3237
  * [12] D. J. Webb and J. van den Berg, “Kinodynamic RRT*: Optimal motion planning for systems with linear differential constraints,” in Proc. IEEE Conf. on Robotics and Automation, 2013.
  * [13] L. Wang, L. Yong and M. H. Ang, ”Hybrid of global path planning and local navigation implemented on a mobile robot in indoor environment,” Proceedings of the IEEE Internatinal Symposium on Intelligent Control, 2002, pp. 821-826.
  * [14] I. Șucan, M. Moll, L. Kavraki, The Open Motion Planning Library, IEEE Robotics & Automation Magazine, 19(4):72–82, December 2012.
  * [15] N. Limpert, S. Schiffer and A. Ferrein. ”A local planner for Ackermann-driven vehicles in ROS SBPL”. 2015 Pattern Recognition Association of South Africa and Robotics and Mechatronics International Conference (PRASA-RobMech), 2015, pp. 172-177.
  * [16] D. Du, D. Hershberger, W. Smart. ”Layered costmaps for context-sensitive navigation”. 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2014.
  * [17] M. Likhachev, G. Gordon, S. Thrun. ”ARA*: Formal Analysis”. School of Computer Science Carnegie Mellon University, 2003.
  * [18] D. Ferguson and M. Likhachev. ”Efficiently using cost maps for planning complex maneuvers”. Lab Papers (GRASP), page 20, 2008.
  * [19] “Generate map with randomly scattered obstacles,” MATLAB. Available: https://www.mathworks.com/help/nav/ref/mapclutter.html. [Accessed: 25-Feb-2022].
  * [20] S. Macenski and I. Jambrecic. ”SLAM Toolbox: SLAM for the dynamic world,” in Journal of Open Source Software, 6(61), 2783.


# Euclidean Invariant Recognition of 2D Shapes Using Histograms of Magnitudes
of Local Fourier-Mellin Descriptors

Xinhua Zhang and Lance R. Williams
Department of Computer Science
University of New Mexico

###### Abstract

Because the magnitude of inner products with its basis functions are invariant
to rotation and scale change, the Fourier-Mellin transform has long been used
as a component in Euclidean invariant 2D shape recognition systems. Yet
Fourier-Mellin transform magnitudes are only invariant to rotation and scale
changes about a known center point, and full Euclidean invariant shape
recognition is not possible except when this center point can be consistently
and accurately identified. In this paper, we describe a system where a
Fourier-Mellin transform is computed at every point in the image. The spatial
support of the Fourier-Mellin basis functions is made local by multiplying
them with a polynomial envelope. Significantly, the magnitudes of convolutions
with these complex filters at isolated points are not (by themselves) used as
features for Euclidean invariant shape recognition because reliable
discrimination would require filters with spatial support large enough to
fully encompass the shapes. Instead, we rely on the fact that normalized
histograms of magnitudes are fully Euclidean invariant. We demonstrate a
system based on the VLAD machine learning method that performs Euclidean
invariant recognition of 2D shapes and requires an order of magnitude less
training data than comparable methods based on convolutional neural networks.

##  1 Introduction

The problem of building a system for 2D shape recognition that is invariant to
rotation, translation and scale change has a long history in computer vision.
Existing approaches to Euclidean invariant recognition have primarily fallen
into three categories: 1) learning invariants by using augmented training
datasets; 2) pooling of non-invariant filter responses; and 3) design of
mathematical invariants. The idea underlying the first category of approaches
is to allow a machine learning system to discover natural invariants by
training it on a dataset of images of different objects augmented with images
of the same objects transformed by rotation, translation and scale change. For
this to work, the augmented training dataset must span the space of different
objects and completely and densely sample the space of Euclidean
transformations. Several examples of work in this category use convolutional
neural networks (CNN) as the machine learning method because these
automatically provide a degree of translation invariance [17, 1, 6], which
lessens the amount of training data required. However, the translation
invariance of a CNN is limited by the lack of translation invariance in its
(non-convolutional) fully-connected layers. Because translation invariance in
practice depends in large part on local pooling operations in these fully-
connected layers, training on a dataset augmented by translations is still
required. In summary, despite the potential of these approaches for very high
recognition rates, the actual degree of Euclidean invariance achieved depends
almost entirely on the number of images in the augmented dataset; denser
sampling of the space of Euclidean transformations yields better invariants
but is often impractical.

Unlike approaches in the first category (which augment the training data in
order to learn invariants) approaches in the second category (in effect)
augment a set of filters. More specifically, they pool the responses of a bank
of non-invariant filters identical up to Euclidean transformation with the
goal of producing a Euclidean invariant net response [42, 38, 20]. For
example, rotating a set of oriented bandpass filters and then averaging the
resulting feature maps yields approximate rotation invariance [38]. Similarly,
applying a filter to an image at different scales yields approximate scale
invariance [20]. As with the data augmentation approaches, the degree of
invariance actually achieved depends on the density and completeness of the
sampling of the space of Euclidean transformations.

As for approaches in the third category, many invariant features are based on
harmonic gratings. This is because the magnitudes of inner products of images
and harmonic gratings [39, 30, 29, 44] are translation invariant. The problem
of designing features invariant to the other components of Euclidean
transformation employ coordinate changes that reduce the transformation to
translation. Rotation invariance, for instance, is translation invariance on
the angular coordinate in polar coordinates [39, 16, 35, 29, 44]. Rotation and
scale invariance can both be reduced to translation invariance in log-polar
coordinates [11, 12, 4, 2, 24, 23, 33]. Interestingly, neuroscience
experiments have shown that neurons in areas V2 and V4 of macaque visual
cortex are selective for harmonic gratings in log-polar coordinates [9, 10,
13]. See Fig. 1. Similar polar separable harmonic gratings can be learned by
machine learning systems trained on sets of images related by a geometric
transformation; when the images in the sets are related by rotation, the
system learns polar harmonic gratings [32]; when the transformations include
both rotation and scaling, the system learns log-polar harmonic gratings [31].

A wide variety of polar separable harmonic grating designs have been used in
systems aimed at Euclidean invariant recognition of 2D shapes. All of these
systems use a harmonic signal in the angular coordinate $\theta$. Since the
angular dimension is periodic, this is the obvious choice. However, the radial
coordinate $r$ is not similarly constrained, and this is where the variety in
the different designs resides. First, the radial coordinate can be logarithmic
or non-logarithmic. The decision to use a logarithmic radial coordinate $\log
r$ is the correct one if scale invariance is required. Second, the function of
the radial coordinate might (or might not) include a harmonic signal.
Including a harmonic signal is useful since it increases the discriminative
ability of the resulting descriptors without compromising scale invariance.
Finally, because the radial dimension is non-periodic, the function of the
radial coordinate should include an envelope to increase locality, and the
function used for this purpose has varied widely. The specific combination of
logarithmic radial variable, harmonic signals in angular and radial variables,
and $\frac{1}{r}$ envelope is embodied in the Fourier-Mellin transform (FMT).
While the $\frac{1}{r}$ envelope provides a weak form of locality while also
preserving scale invariance, it doesn’t satisfy the stronger locality
constraint defined by [30]. Conversely, while a polynomial envelope satisfies
the stronger locality constraint, it would seem to be incompatible with scale
invariance. Significantly, in this paper we show that log-polar harmonic
filters with polynomial envelopes can (in fact) be used for scale invariant
recognition if the vector of responses is appropriately normalized.

There are other approaches to achieving translation invariance distinct from
harmonic gratings. Methods based on histograms (or other statistics) of
convolutions are inherently translation invariant because spatial information
is discarded. Representative methods assume that inputs have a multimodal
distribution and computes a visual vocabulary to model this distribution. For
example, the Fisher vector method [34] uses Gaussian mixture models as its
visual vocabulary, while the bag of words [40] and vector of locally
aggregated descriptor (VLAD) [19] methods use centers computed by k-means
clustering.

In this paper we propose to use the statistics of locally rotation and scale
invariant features to construct globally Euclidean invariant descriptors.
Using feature statistics, e.g., of wavelet subbands, rather than the features
themselves, is an approach more commonly used in texture analysis. We extend
this idea and show that it can also be used for Euclidean invariant 2D shape
recognition. Specifically, we apply the VLAD machine learning algorithm to the
magnitudes of local Fourier-Mellin descriptors—log-polar separable harmonic
gratings with polynomial envelopes. We explore the quality of locality for
different choices of envelope and show that descriptors based on polynomial
envelopes exhibit superior scale invariance and result in better image
reconstructions from magnitudes alone.

Figure 1: Log-polar harmonic gratings. The phases of complex values are mapped
to hue with red indicating positive real. Each row is a distinct radial
frequency. From top to bottom: 0, 1, 2. Each column is a distinct angular
frequency. From left to right: -5, -3, -1, 0, 1, 3, 5.

##  2 Related Work

Chi et al. [2] use signatures based on the energy of orthogonal wavelet
transforms of log-polar images as Euclidean invariants. The technique they
describe uses an adaptive row shift-invariant wavelet transform to achieve a
degree of translation invariance in the log-polar domain because (unlike
harmonic gratings) discrete wavelet transforms are not intrinsically
translation invariant.

Esteves et al. [6] train a CNN to recognize objects in log-polar transformed
images. However, a separate CNN is used to learn the object center points
before the log-polar transform is applied. Finally, as with all CNNs, the
translation invariance achievable without extra training data is limited, and
the rotation invariant recognition rate would be low if the system were not
trained on augmented datasets.

Worrall et al. [44] uses convolutions with harmonic gratings in polar
coordinates as input to a CNN to achieve rotation and translation invariant
recognition of 2D shapes. However, they did not address the question of scale
invariance and rotation and translation invariance are only achieved by
training on augmented datasets. Pew et al. [35] use a set of orthogonal basis
functions in polar coordinates as rotation invariant features for 2D shape
recognition. Like [44], they do not address the issue of scale invariance.

The prior work most similar to our own is by Gotze et al. [12], who used
rotation and translation invariant descriptors based on magnitudes of inner
products with polar separable harmonic gratings with Gaussian envelopes. To
avoid the problem of identification of the center point, these descriptors are
computed at every point in the image using convolution. The dimensionality of
the descriptor computed at each point is relatively low, and unlike the system
we describe, there is no use of the statistics of descriptors within extended
regions. Consequently, acceptable recognition accuracy can only be achieved by
making the spatial support of the filters large enough to fully encompass the
shape being recognized. Finally, unlike our system, which uses polynomial
envelopes to achieve scale invariance, the use of the Gaussian envelope by
[12] makes this impossible.

Lehiani et al. [27] describe a system that uses magnitudes of analytical
Fourier-Mellin transform computed at keypoints as a Euclidean invariant
signature. Like [12], the dimensionality of the descriptor at each keypoint is
relatively low and acceptable recognition accuracy is only achieved by using
Fourier-Mellin basis functions with large spatial support.

Hoang et al. [14, 15] describe a Euclidean invariant descriptor which combines
the radon transform and magnitudes of a 1D Fourier-Mellin transform. While the
basic idea of constructing an invariant descriptor by converting rotation and
scaling to translation is similar to ours, their descriptor is not spatially
localized and our results demonstrate the importance of localization.

Finally, Mellor, Hong and Brady [30] used histograms of responses of polar
separable filters with a polynomial envelope for texture recognition. Our work
is different from theirs in three ways: 1) The filters they use are not
harmonic gratings and are unrelated to Fourier-Mellin basis functions; 2) They
use the joint histogram of only two Euclidean invariant features while we use
VLAD (a non-parametric generalization of Fisher vectors) to represent the
joint statistics of a large number of Euclidean invariant features; 3) They
only applied their method to Euclidean invariant recognition of textures,
while our method is also applied to Euclidean invariant recognition of 2D
shapes.

##  3 Euclidean Invariant Descriptors

An output representation is equivariant with respect to a geometric
transformation if the operation that produces it commutes with the
transformation. For example, convolution commutes with rotation and
translation of the input image if the filter is isotropic. In contrast, an
output representation is invariant with respect to a geometric transformation
if it is unchanged by it. For example, histograms are unchanged by rotation
and translation of the input image.

Figure 2: An illustration of the relationship between local transforms and
global transforms. When a patch is rotated and scaled with with respect to the
origin, all points of the patch are rotated and scaled by the same amount with
respect to the center of the patch.

###  3.1 From Local to Global Invariance

Consider an image $I$ with a gray value histogram $H$. A Euclidean
transformation ${\cal T}$ of $I$ about its origin consists of a scale change
$s$, rotation by an angle $\phi$, and translation by a vector, $[u,v]^{\rm
T}$. This transformation maps the image and its histogram as follows:

$\displaystyle I\left(\left[\begin{array}[]{c}x\\\ y\end{array}\right]\right)$
$\displaystyle\stackrel{{\scriptstyle\cal T}}{{\rightarrow}}$ $\displaystyle
I\left(\left[\begin{array}[]{c}s\>(\>x\>\cos\phi+\>y\>\sin\phi)+u\\\
s\>(\>y\>\cos\phi-\>x\>\cos\phi)+v\end{array}\right]\right)$ $\displaystyle
H(i)$ $\displaystyle\stackrel{{\scriptstyle\cal T}}{{\rightarrow}}$
$\displaystyle s^{2}\;H(i).$

Because scaling an image by $s$ increases the size of areas by $s^{2}$, and
because the histogram represents the area possessing a single gray value, $H$
is multiplied by $s^{2}$. It follows that $H\>/\>||H||_{2}$ is invariant with
respect to Euclidean transformation.

It is possible to construct rotation and scale invariant features with spatial
support greater than one pixel. A Fourier-Mellin basis function
$f_{\omega_{r},\omega_{\theta}}$ is a separable function in log-polar
coordinates:

$f_{\omega_{r},\omega_{\theta}}=e^{-j\omega_{r}\log
r}e^{-j\omega_{\theta}\theta}r^{-1}.$

The magnitude of the inner product of a Fourier-Mellin basis function and a
circular image patch $p$ is invariant to rotation and scaling of the patch
because these are translation in log-polar coordinates and magnitude is not
changed by translation:

$|\left<f_{\omega_{r},\omega_{\theta}},\>p(\log
r,\theta)\right>|=|\left<f_{\omega_{r},\omega_{\theta}},\>p(\log r\\!+\\!\log
s,\theta\\!+\\!\phi)\right>|.$

In an image consisting of a set $P$ of overlapping circular patches, a patch
$p$ is accessed by specifying its global position in Cartesian coordinates
$[x,y]^{\rm T}$. Gray values are accessed by specifying their local positions
within $p$ in log-polar coordinates $[\log r,\theta]^{\rm T}$. Now, consider
the effect of the Euclidean transformation ${\cal T}$ on the set of patches
and a histogram $H_{\omega_{r},\omega_{\theta}}$ of the magnitudes of inner
products of the patches with a Fourier-Mellin basis function
$f_{\omega_{r},\omega_{\theta}}$:

$\displaystyle P\left(\left[\begin{array}[]{c}x\\\ y\\\ \hline\cr\log r\\\
\theta\end{array}\right]\right)$ $\displaystyle\stackrel{{\scriptstyle\cal
T}}{{\rightarrow}}$ $\displaystyle
P\left(\left[\begin{array}[]{c}s\>(\>x\>\cos\phi+y\>\sin\phi)+u\\\
s\>(\>y\>\cos\phi-x\>\cos\phi)+v\\\ \hline\cr\log r+\log s\\\
\theta+\phi\end{array}\right]\right)$ $\displaystyle
H_{\omega_{r},\omega_{\theta}}(i)$ $\displaystyle\stackrel{{\scriptstyle\cal
T}}{{\rightarrow}}$ $\displaystyle s^{2}\;H_{\omega_{r},\omega_{\theta}}(i).$

Although the patch centers move, the magnitudes of inner products with
$f_{\omega_{r},\omega_{\theta}}$ are invariant with respect to the local
geometric transformation induced by ${\cal T}$ for all patches because this is
just translation by $[\log s,\phi]^{\rm T}$. Consequently, magnitude of
convolution with $f_{\omega_{r},\omega_{\theta}}$ commutes with the global
transformation ${\cal T}$. It follows that
$|I*f_{\omega_{r},\omega_{\theta}}|$ is equivariant with respect to Euclidean
transformation. Finally, because $|I*f_{\omega_{r},\omega_{\theta}}|$ is
equivariant with respect to Euclidean transformation,
$H_{\omega_{r},\omega_{\theta}}\>/\>||H_{\omega_{r},\omega_{\theta}}||_{2}$ is
invariant with respect to Euclidean transformation. See Fig. 2.

###  3.2 Local Fourier-Mellin Transform

Given an image $I(r,\theta)$ in polar coordinates, the local Fourier-Mellin
transform $LFMT(\omega_{r},\omega_{\theta},\alpha)$ is defined as

$\frac{1}{2\pi}\int_{0}^{2\pi}\int_{0}^{R}I(r,\theta)e^{-j\omega_{r}\log
r}e^{-j\omega_{\theta}\theta}r^{\alpha}drd\theta$

where $\omega_{r}$ and $\omega_{\theta}$ are radial and angular frequencies,
$R$ is the range of the transform, and $\alpha$ defines the envelope shape.
Fourier-Mellin transform (FMT) is the case when $\alpha=-1$. The envelope is
scale invariant because the Mellin transform implies
$\frac{dr}{r}=\frac{d(sr)}{sr}$. Therefore, the magnitude of FMT is scale
invariant:

$\begin{split}&\quad\left|\>\frac{1}{2\pi}{\textstyle\int_{0}^{2\pi}\int_{0}^{sR}I\left(r/s,\theta\right)e^{-j\omega_{r}\log
r}e^{-j\omega_{\theta}\theta}\frac{dr}{r}d\theta}\>\right|\\\
&=\left|\>\frac{1}{2\pi}{\textstyle\int_{0}^{2\pi}\int_{0}^{R}I\left(r^{\prime},\theta\right)e^{-j\omega_{r}\log
r^{\prime}}e^{-j\omega_{\theta}\theta}\frac{dr^{\prime}}{r^{\prime}}d\theta}\>\right|\end{split}$

where $r^{\prime}=\frac{r}{s}$ and $|\cdot|$ is magnitude. For other values of
$\alpha$, there is an extra term $s^{\alpha+1}$, because of
$(sr)^{\alpha}d(sr)$. Consequently, the magnitude of the transform is not
scale invariant. Yet the L2 normalized vector comprised of magnitudes at all
frequencies is. Several works [39, 30, 29, 44] use relative phases instead of
magnitudes because relative phases are both rotation and contrast invariant.
We use magnitudes in our work because the number of relative phases is
$M^{2}\times N^{2}$, and this is quite large. Since the vector of magnitudes
can be made contrast invariant by dividing it by its L2 norm, magnitudes are
the better choice.

In order to compute the LFMD, we need to derive a set of convolution kernels
by discretely sampling log-polar harmonic gratings with polynomial envelopes.
Several issues need to be addressed. First, the kernels should have odd size
so that the center of the kernel is pixel-centered and because $\log
0=-\infty$, the center of the kernel should be 0. Second, because the kernels
are constructed by discrete sampling of log-polar harmonic gratings, the
problem of aliasing in the $\theta$ dimension becomes significant as the
angular frequency $\omega_{\theta}$ increases. The solution would seem to be
to modify the envelope function by zeroing the kernel within a radius large
enough to prevent aliasing of a given angular frequency.111 Interestingly,
CNNs trained on rotation augmented datasets seem to discover the same
solution, i.e., functions resembling polar harmonic gratings with holes in
their centers [32, 31]. Accordingly, we define a log-polar harmonic grating
with polynomial envelope and a hole in the center as follows

${\small h_{\omega_{r},\omega_{\theta},\alpha}(x,y)=}$
${\small\left\\{\begin{array}[]{ll}0&\\!\\!{\rm
if}\sqrt{x^{2}+y^{2}}\leq\sigma\omega_{\theta}\\\
\frac{1}{2\pi}(\sqrt{x^{2}+y^{2}})^{(\alpha-j\omega_{r})}e^{-j\omega_{\theta}{\rm
atan}\left(\frac{x}{y}\right)}&\\!\\!\mathrm{otherwise.}\end{array}\right.}$

where the $\exp(-j\omega_{r}\log r)$ has been replaced by $r^{-j\omega_{r}}$.
Experiments show that $\sigma=\frac{2}{\pi}$ gives the best performance.
Finally, given an image $I$, a local Fourier-Mellin descriptor $LFMD_{\alpha}$
at pixel $(i,j)$ is defined as

$n\left(\left[|\\{I*h_{\omega_{r_{1}}\\!,\omega_{\theta_{1}}\\!,\alpha}\\}(i,j)|,\dots,|\\{I*h_{\omega_{r_{M}}\\!,\omega_{\theta_{N}}\\!,\alpha}\\}(i,j)|\right]^{\rm{T}}\right)$

where $n(\vec{x})=\vec{x}\;/\;\|\vec{x}\|_{2}$.

###  3.3 Envelope Locality and Image Reconstruction

Different envelopes are used for different purposes. The $\alpha=-1$ envelope
of the FMT is best for scale invariance. Analytical Fourier-Mellin transform
(AFMT) [11, 4] is the case when $\alpha>-1$. According to Mellor, Hong and
Brady [30], a filter $g(r)$ is localized if

$\int_{0}^{R}2\pi rg(r)^{2}dr>\int_{R}^{\infty}2\pi rg(r)^{2}dr$

for all $R>0$. Envelopes satisfy this constraint when $\alpha<-1$.
Consequently, the envelopes used in FMT and AFMT are not localized. The effect
of different values of $\alpha$ on the appearance of the filters is shown in
Fig. 3. Because the pros and cons of the choice of envelope shape are unclear,
we evaluate the effect of envelope shape empirically by gauging performance on
two problems: image classification (Section 4) and image reconstruction
(Section 5).

For invertible transforms, such as the polar harmonic transform [35] or the
AFMT [4], the original image can be recovered simply by applying the inverse
transform.222Of course if the basis functions of the transform are localized,
and the image size is larger than the spatial support of these functions, then
it will be unable to recover the entire image. However, there is no way to
recover the original image from the transform coefficient magnitudes alone.
Since these are the actual Euclidean invariants that will be used for 2D shape
recognition, the question of whether or not they are sufficient by themselves
for reconstructing the original image is not purely academic. It has been
previously shown that images can sometimes be reconstructed from the
magnitudes of values resulting by convolution with complex-valued filters
[37]. Although there is a sign ambiguity in the reconstruction, and this is
unavoidable since magnitudes discard sign, the qualities of the
reconstructions are otherwise good. The results described in [37] uses the
magnitudes of complex Gabor wavelets. Because they are the products of
harmonic gratings in Cartesian coordinates and Gaussian envelopes, these
filters possess local translation invariance. In this paper, we show that
images can be reconstructed from the magnitudes of LFMDs, the products of
harmonic gratings in log-polar coordinates and polynomial envelopes, and which
possess local rotation and scale invariance.

Given an image $I$, a reconstructed image $\tilde{I}$, and a local Fourier-
Mellin filter $h_{\omega_{r},\omega_{\theta},\alpha}$, the reconstruction
error $E$ is defined as follows

$\begin{split}\sum_{\omega_{r}}\sum_{\omega_{\theta}}{\left(\left|\frac{I*h_{\omega_{r},\omega_{\theta},\alpha}\phantom{\tilde{I}}}{\||I*h_{\omega_{r},\omega_{\theta},\alpha}|\|_{2}}\right|^{2}-\>\>\left|\frac{\tilde{I}*h_{\omega_{r},\omega_{\theta},\alpha}}{\||\tilde{I}*h_{\omega_{r},\omega_{\theta},\alpha}|\|_{2}}\right|^{2}\right)^{2}}.\end{split}$

At every position, the convolution is divided by the L2 norm of magnitudes
across frequencies. The reconstruction $\tilde{I}$ is initialized randomly and
updated by gradient descent where
$\Delta\tilde{I}(x,y)=\frac{dE}{d\tilde{I}(x,y)}$. The normalization is
necessary for LFMDs to be scale and contrast invariant. Although the
definition of error in [37] does not have the normalization term, it is
otherwise the same as the one we used in the present work.

(a)

(b)

(c)

(d)

Figure 3: Effects of $\alpha$ on Local Fourier-Mellin filters. Filter size is
$64\times 64$. (a) $\alpha$ = -0.5 (b) -1 (c) -2 (d) -3.

###  3.4 LFMD-VLAD

Besides standard VLAD [19], we also use the power-law normalization. More
specifically, given a vector $\vec{v}=[v_{1},\dots,v_{N}]^{\rm T}$, a
normalized vector $\vec{v}^{\prime}$ equals
$[v^{\prime}_{1},\dots,v^{\prime}_{N}]^{\rm T}$ where
$\vec{v}^{\prime}_{i}=|v_{i}|^{\beta}\times\mathrm{sign}(v_{i})$ [18]. The
scale, translation and contrast invariant descriptor, LFMD-VLAD, is computed
using the following procedure:

  1. 1.

Compute a $MN$ dimensional LFMD at every point.

  2. 2.

Compute $K$ centers, $\vec{c}_{k}$, of LFMDs using $k$-means clustering.

  3. 3.

Compute
$\vec{v}_{k}\\!=\\!{\textstyle\sum_{\\{\vec{x}\>|\>\vec{n}(\vec{x})=\vec{c}_{k}\\}}{\vec{x}-\vec{c}_{k}}}$
where $\vec{n}(\vec{x})$ is the nearest neighbor of $\vec{x}$.

  4. 4.

Construct a $K\times MN$ dimensional vector $\vec{V}$ by concatenating the
$\vec{v}_{k}$.

  5. 5.

Apply $\beta=0.5$ power law normalization followed by L2 normalization to
$\vec{V}$.

###  3.5 Scale Invariance and Object Background Ratio

When an image of an object is scaled, the ratio of the area of the object and
the area of the background changes. This affects the LFMD magnitude histogram:

$\displaystyle H(0)$ $\displaystyle\rightarrow$ $\displaystyle
H(0)+{\textstyle\sum_{i>0}\left[H(i)-s^{2}H(i)\right]}$ $\displaystyle H(i)$
$\displaystyle\rightarrow$ $\displaystyle s^{2}H(i).$

where $H$ is the histogram of an image of an isolated object against a
uniformly zero valued background. In a case such as the above, the problem can
be solved by thresholding. However, the problem is more complex in natural
images, and other techniques (beyond the scope of the present work) are
needed.

##  4 Image Classification

In this section, we gauge the degree to which the LFMD-VLAD system achieves
Euclidean invariant recognition by testing it in three different problem
domains: 1) the MNIST database [26] for handwritten digit recognition; 2) the
Swedish leaves database [41] for 2D shape recognition; and 3) the KTH-TIPS
database [8] for texture recognition.

We apply rotation and scaling to the MNIST and Swedish leaves databases. All
transforms use bicubic interpolation. Our classifier is trained on original
images and tested on transformed images. We did not augment the set of
training images in anyway, e.g., by inclusion of additional images related to
the training set by translation, rotation and scaling. We demonstrate that
LFMD-VLAD does not need to be trained on augmented datasets to achieve
Euclidean invariant recognition. In all classification experiments,
$\omega_{r}\in\\{0,1,2\\}$ and $\omega_{\theta}\in\\{-5,-4,\dots,5\\}$. The
dimensionality of the descriptor, $MN$, is therefore $3\times 11=33$. There
are $K=64$ centers computed by the $k$-means clustering algorithm so the
dimension of the LFMD-VLAD system is $MN\times K=33\times 64=2112$. For color
images, the same procedure is applied to the three channels independently, so
the descriptor dimensionality is three times that of gray scale images. We use
the linear SVM described in [7] with the parameter $C=1$. Several different
values of $\alpha$ are evaluated in order to ascertain the effect of using
different envelopes. When $\alpha$ equals 0, the filters are non-localized
log-polar harmonic gratings; when $\alpha$ equals 0.5, the filters implement
the analytical FMT described by [4]; when $\alpha$ equal to -1, the filters
are the basis functions of the classical FMT. Finally, we will show
experimentally that setting $\alpha$ equal to -2 or -3 gives the best locality
without sacrificing scale invariance.

We use the algorithm described by [6] as a benchmark to compare with our own
on the MNIST database [26]. Given an image, it is first log-polar transformed
into a $256\times 256$ image using the center of the image. Afterwards, an
AlexNet [25] with $224\times 224$ random cropping and random flipping is
trained on the log-polar transformed images. The central idea of this approach
is to convert rotation and scaling in Cartesian into translation in log-polar,
and CNNs can build up translation-invariance by local pooling layer after
layer. Regarding the methods proposed in [6], we do not use the same CNNs
because they are not better than AlexNet with respect to translation-
invariance in theory. We do not use the origin predictor either, because our
objective is not to maximize recognition accuracy, but to illustrate the
degree to which the CNN depends on the use of augmented training sets to
achieve Euclidean invariance.

Table 1: MNIST classification error rates ($\%$) using a training set which consists of only original images. For testing sets, O is original, R is rotated, S is scaled and R+S is rotated and scaled. Names for methods and training sets are in the first column. For training sets, O is original, L is log-polar transformed and + indicates an R+S augmented training set. | $\alpha$ | O | R | S | R + S
---|---|---|---|---|---
AlexNet (O) | – | 1.27 | 52.31 | 7.33 | 58.48
AlexNet (O+) | – | 5.28 | 5.31 | 6.04 | 6.14
AlexNet (L) | – | 1.09 | 37.56 | 2.03 | 41.44
AlexNet (L+) | – | 3.87 | 3.75 | 3.91 | 3.72
LFMD-VLAD | 0 | 5.01 | 5.01 | 5.08 | 5.03
LFMD-VLAD | -0.5 | 3.92 | 3.9 | 4.12 | 4.01
LFMD-VLAD | -1 | 2.75 | 2.78 | 3.09 | 3.12
LFMD-VLAD | -2 | 2.46 | 3.02 | 5.98 | 6.7
LFMD-VLAD | -3 | 7.04 | 13.34 | 21.26 | 27.11

###  4.1 MNIST

The standard MNIST database contains 60,000 gray scale images of digits of
size $28\times 28$.333We omitted the digit ‘9’ since the purpose of our
research is Euclidean invariant recognition and the ‘9’ cannot be
distinguished from ‘6’ when rotation is permitted. To create the augmented
datasets used for testing (and training also in the case of the CNN), each
image was resized to $48\times 48$, geometrically transformed in one of the
following ways:

  1. 1.

Original: No transformation.

  2. 2.

Rotation: A random rotation between $0$ and $2\pi$

  3. 3.

Scaling: A random scaling between 0.5 and 1.5

  4. 4.

Rotation and scaling: The combination of the rotation and scaling
transformations defined above.

After the geometric transformation was applied, the image was padded to size
$96\times 96$. In this way, four sets of images were created. We do not
augment the dataset in like [6] because we want to show that the proposed
descriptor can handle transformations in a wider range. Also, this can
decrease performance for CNNs more than using a smaller range. For every
transform, the training set has 9,000 images (1,000 per category) images, and
the testing set has 45,000 images (5,000 per category) images. Descriptors
used by $k$-means are extracted from the first 1,000 images with stride 4, but
those used for constructing LFMD-VLAD for one image are extracted with stride
1.

The results are shown in Table 1. AlexNet fails when used with log-polar
filters when trained on a non-augmented dataset because the CNN does not
provide enough translation invariance. The fact that it is more scale
invariant than rotation invariant suggests that it is unable to learn the
periodic nature of the $\theta$ dimension. Although it is able to learn
rotation and scale invariance when trained on an augmented dataset, this
requires approx. 10 times more training data than is required by the LFMD-VLAD
approach (Fig 4).

Figure 4: Effects of training set size. The format of the legend is algorithm,
training set, testing set, where O is original and R+S is rotated and scaled.
The dash line marks the best error rate achieved by the CNN when trained on
the augmented dataset.

###  4.2 Swedish leaves

The Swedish leaves database consists of 15 categories and 75 images per
category and is primarily used for 2D shape recognition. In our experiments,
375 images (25 per category) are randomly chosen to serve as the training
dataset and the remaining 750 images serve as the testing dataset. The same
geometric transforms that were used in the MNIST dataset experiment were
applied to construct the augmented testing dataset except that the scale
factor is between 0.5 and 1. Afterwards, images are padded to size $256\times
256$. The process of creating training and testing datasets is repeated five
times and the reported classification accuracies are the average values
achieved. The descriptors used for $k$-means clustering are extracted from all
training images with stride 4.

Although the original leaf images appear to have a white background, the
background is not uniform. To remedy this, the images are preprocessed as
follows: 1) gray values are rescaled to a range between 0 and 255 and then
inverted; 2) gray values less than the average for the image are set to zero.
This results in leaf images with uniformly black backgrounds.

Data augmentation sometimes makes AlexNet’s performance worse. Since the R+S
augmentation only transforms every single image once, it always reduces the
number of similar images between training set and testing set. The fact that
there are only 375 training images makes this phenomenon more obvious than in
the MNIST experiment. We observe that the optimum value of $\alpha$ is
different than it was for the MNIST experiment, with more localized filters
($\alpha=-3$) yielding increased accuracy. Both shape recognition algorithms
[21, 36] use rotation-invariance to facilitate shape recognition. However,
neither shows how much invariance is actually achieved and LFMD-VLAD’s
recognition accuracy is competitive with both on the leaf dataset.

Table 2: Swedish leaves classification accuracy ($\%$) using a training set which consists of only original images. For testing sets, O is original, R is rotated, S is scaled and R+S is rotated and scaled. Names for methods and training sets are in the first column. For training sets, O is original, L is log-polar transformed and + indicates an R+S augmented training set. | $\alpha$ | O | R | S | R + S
---|---|---|---|---|---
[21] | – | 95.33 | – | – | –
[36] | – | 99.38 | – | – | –
AlexNet (O) | – | 98.13 | 28.40 | 75.67 | 24.72
AlexNet (O+) | – | 92.51 | 40.75 | 89.76 | 39.25
AlexNet (L) | – | 75.15 | 74.75 | 79.47 | 82.37
AlexNet (L+) | – | 73.76 | 78.21 | 80.11 | 81.83
LFMD-VLAD | 0 | 91.46 | 84.66 | 79.90 | 78.2
LFMD-VLAD | -0.5 | 92.8 | 85.84 | 82.04 | 80.8
LFMD-VLAD | -1 | 94.18 | 88.6 | 84.28 | 82.52
LFMD-VLAD | -2 | 98.14 | 96.6 | 78.96 | 78.48
LFMD-VLAD | -3 | 99.14 | 98.40 | 95.92 | 95.84

###  4.3 KTH-TIPS

KTH-TIPS contains 10 categories and 81 texture images per category. Images
from one category vary in scale, orientation and contrast. For this database,
we do not need to apply any geometric transformations because they are already
included in the dataset. For cross validation purposes, training datasets are
created in different sizes: 50 (5 per category), 200 (20 per category) and 400
(40 per category). The remaining images are used to define the testing
dataset. Images are chosen at random and this process is repeated five times.
The classification accuracy is the average of these five repetitions. Because
textures have stationary spatial statistics, the issue with uniform background
is not relevant for this dataset. Consequently, no thresholding of the LFMD
values was used.

The algorithm proposed in [38] is state-of-the-art for the KTHTIPS dataset. It
outperforms LFMD-VLAD in the 20 per category and 40 per category cases.
However, LFMD-VLAD has the best performance in the 5 per category case.
Increasing the number of training images does not improve LFMD-VLAD’s
performance as much it does in the $\alpha=-2$ or $\alpha=-3$ cases. This
suggests a weak relationship between Euclidean-invariance and texture
recognition rate.

Table 3: KTHTIPS classification accuracy ($\%$). Names for methods and training sets are in the first column. O stands for trained on original images and L stands for trained on log-polar transformed images.  | $\alpha$ | 5 | 20 | 40
---|---|---|---|---
[38] | – | 84.3 | 98.3 | 99.4
AlexNet (O) | – | 79.29 | 86.98 | 93.42
AlexNet (L) | – | 57.92 | 73.02 | 78.83
LFMD-VLAD | 0 | 66.1 | 87.8 | 88.8
LFMD-VLAD | -0.5 | 66.7 | 86.0 | 89.3
LFMD-VLAD | -1 | 71.5 | 88.3 | 90.3
LFMD-VLAD | -2 | 85.5 | 95.2 | 95.4
LFMD-VLAD | -3 | 85.3 | 94.8 | 95.7

##  5 Image Reconstruction

In this section, we demonstrate the reconstruction of an image using LFMDs and
the algorithm described in Sec. 3.3. We use the same filters as in Sec. 4,
except that there is no hole at the center of the filter. In Eq. (3.2), the
value is set to 0 if $\sqrt{x^{2}+y^{2}}\leq\sigma\omega_{\theta}$. For the
sake of reconstruction quality, it is set to 0 only if $\sqrt{x^{2}+y^{2}}=0$.
Although a hole at the center can reduce aliasing and improve classification
accuracy, the missing information at the center results in noisy
reconstructions. Consequently, to better illustrate the idea that images can
be reconstructed from LFMDs, we use descriptors with slightly decreased scale
invariance. Finally, we observe that the phase of the reconstruction depends
on the initial value image and this image is random. See Fig. 5(c), 5(d) and
Fig. 5(e), 5(f).

The FMT fails to reconstruct the image (Fig. 5(b)). In fact, reconstruction
always fails when $\alpha\geq-1$. This demonstrates that the locality of the
filters is critical to reconstruction. Best reconstruction quality (Fig. 5(c))
is achieved when $\alpha=-2$. However, because the filter support is very
small when $\alpha=-3$, it is difficult to achieve consistency in the
reconstructed local phases (Fig. 5(f)). Nevertheless, it is clear that
translation invariant descriptors, such as magnitudes of Gabor wavelets [37],
and rotation and scale invariant descriptor, such as LFMDs, are both able to
reconstruct images. The phase information (presumed lost) is implicit in the
magnitudes of the responses of the complex filters if the filters are
sufficiently well localized.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: (a) Original image of size $256\times 256$. (b) Poor reconstruction
for $\alpha=-1$. (c) Improved reconstruction for $\alpha=-2$. (d) Alternative
reconstruction for $\alpha=-2$. (e) Poor reconstruction for $\alpha=-3$. (f)
Alternative for $\alpha=-3$.

##  6 Conclusion

This paper introduced a locally rotation and scale invariant descriptor called
LFMD. The rotation and scale invariance of LFMDs computed for local patches
results in rotation, translation and scale equivariance of non-local image-
based representations comprised of them. This equivariance, in turn, results
in rotation, translation and scale invariant 2D shape recognition when
combined with VLAD. Prior work on 2D shape recognition has ignored scale,
using image-based representations that are solely translation and rotation
equivariant [28, 43, 44, 5, 3, 22].

We showed that a system based on LFMD-VLAD is capable of Euclidean invariant
2D shape recognition without the necessity of training on augmented datasets.
We also showed that images can be reconstructed from LFMDs alone. This is
similar to what others have done using magnitudes of convolutions with Gabor
wavelets. Results on 1) Euclidean invariant 2D shape recognition; and 2)
recovery of phases from magnitudes in images convolved with complex filters
suggest that filter locality is essential.

In future work, we plan to generalize the method described in this paper to
images with cluttered backgrounds and build a hierarchical network exploiting
the Euclidean equivariance of non-local image-based representations based on
LFMDs.

Acknowledgement Xinhua Zhang gratefully acknowledges the support of the New
Mexico Consortium.

## References

  * [1] G. Cheng, P. Zhou, and J. Han.  RIFD-CNN: Rotation-Invariant and Fisher Discriminative Convolutional Neural Networks for Object Detection.  In CVPR, 2016.
  * [2] Chi-Man Pun and Moon-Chuen Lee.  Log-polar wavelet energy signatures for rotation and scale invariant texture classification.  TPAMI, 2003.
  * [3] T. S. Cohen and M. Welling.  Group Equivariant Convolutional Networks.  In ICML, 2016.
  * [4] S. Derrode and F. Ghorbel.  Robust and Efficient Fourier–Mellin Transform Approximations for Gray-Level Image Reconstruction and Complete Invariant Description.  Computer Vision and Image Understanding, 2001.
  * [5] S. Dieleman, J. De Fauw, and K. Kavukcuoglu.  Exploiting Cyclic Symmetry in Convolutional Neural Networks.  In ICML, 2016.
  * [6] C. Esteves, C. Allen-Blanchette, X. Zhou, and K. Daniilidis.  Polar Transformer Networks.  In ICLR, 2018.
  * [7] R.-E. Fan, K.-W. Chang, C.-J. Hsieh, X.-R. Wang, and C.-J. Lin.  LIBLINEAR: A Library for Large Linear Classification.  JMRL, 2008.
  * [8] M. Fritz, E. Hayman, B. Caputo, and J.-o. Eklundh.  The KTH-TIPS database.  2004\.
  * [9] J. Gallant, J. Braun, and D. Van Essen.  Selectivity for polar, hyperbolic, and Cartesian gratings in macaque visual cortex.  Science, 1993.
  * [10] J. L. Gallant, C. E. Connor, S. Rakshit, J. W. Lewis, and D. C. Van Essen.  Neural responses to polar, hyperbolic, and Cartesian gratings in area V4 of the macaque monkey.  Journal of Neurophysiology, 1996.
  * [11] F. Ghorbel.  A complete invariant description for gray-level images by the harmonic analysis approach.  Pattern Recognition Letters, 1994.
  * [12] N. Gotze, S. Drue, and G. Hartmann.  Invariant object recognition with discriminant features based on local fast-Fourier Mellin transform.  In ICPR, 2000.
  * [13] J. Hegdé and D. C. Van Essen.  Selectivity for complex shapes in primate visual area V2.  The Journal of Neuroscience, 2000.
  * [14] T. V. Hoang and S. Tabbone.  A geometric invariant shape descriptor based on the radon, fourier, and mellin transforms.  ICPR, 2010.
  * [15] T. V. Hoang and S. Tabbone.  Invariant pattern recognition using the RFM descriptor.  Pattern Recognition, 2012.
  * [16] G. Jacovitti and A. Neri.  Multiresolution circular harmonic decomposition.  IEEE Transactions on Signal Processing, 2000.
  * [17] M. Jaderberg, K. Simonyan, A. Zisserman, and K. Kavukcuoglu.  Spatial Transformer Networks.  In NIPS, 2015.
  * [18] H. Jégou and O. Chum.  Negative evidences and co-occurences in image retrieval: The benefit of PCA and whitening.  In ECCV, 2012.
  * [19] H. Jegou, M. Douze, C. Schmid, and P. Perez.  Aggregating local descriptors into a compact image representation.  In CVPR, 2010.
  * [20] A. Kanazawa, A. Sharma, and D. Jacobs.  Locally Scale-Invariant Convolutional Neural Networks.  ArXiv e-prints, 2014.
  * [21] Q. Ke and Y. Li.  Is Rotation a Nuisance in Shape Recognition?  In CVPR, 2014.
  * [22] J. J. Kivinen and C. K. Williams.  Transformation equivariant Boltzmann machines.  In ICANN, 2011.
  * [23] I. Kokkinos, M. M. Bronstein, and A. Yuille.  Dense Scale Invariant Descriptors for Images and Surfaces.  Research Report, 2012.
  * [24] I. Kokkinos and A. Yuille.  Scale invariance without scale selection.  In CVPR, 2008.
  * [25] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  ImageNet Classification with Deep Convolutional Neural Networks.  In NIPS, 2012.
  * [26] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 1998.
  * [27] Y. Lehiani, M. Maidi, M. Preda, and F. Ghorbel.  Image invariant description based on local Fourier-Mellin transform.  In ICSIPA, 2017.
  * [28] K. Lenc and A. Vedaldi.  Understanding Image Representations by Measuring Their Equivariance and Equivalence.  IJCV, 2018.
  * [29] K. Liu, H. Skibbe, T. Schmidt, T. Blein, K. Palme, T. Brox, and O. Ronneberger.  Rotation-Invariant HOG Descriptors Using Fourier Analysis in Polar and Spherical Coordinates.  IJCV, 2014.
  * [30] M. Mellor, Byung-Woo Hong, and M. Brady.  Locally Rotation, Contrast, and Scale Invariant Descriptors for Texture Analysis.  TPAMI, 2008.
  * [31] R. Memisevic.  Gradient-based learning of higher-order image features.  In ICCV, 2011.
  * [32] R. Memisevic and G. E. Hinton.  Learning to Represent Spatial Transformations with Factored Higher-Order Boltzmann Machines.  Neural Computation, 2010.
  * [33] J. Mennesson, C. Saint-Jean, and L. Mascarilla.  Color Fourier–Mellin descriptors for image recognition.  Pattern Recognition Letters, 2014.
  * [34] F. Perronnin and C. Dance.  Fisher Kenrels on Visual Vocabularies for Image Categorizaton.  CVPR, 2006.
  * [35] Pew-Thian Yap, Xudong Jiang, and A. C. Kot.  Two-Dimensional Polar Harmonic Transforms for Invariant Image Representation.  TPAMI, 2010.
  * [36] X. Qi, R. Xiao, C.-G. Li, Y. Qiao, J. Guo, and X. Tang.  Pairwise Rotation Invariant Co-Occurrence Local Binary Pattern.  TPAMI, 2014.
  * [37] L. Shams and C. von der Malsburg.  The role of complex cells in object recognition.  Vision Research, 2002.
  * [38] L. Sifre and S. Mallat.  Rotation, Scaling and Deformation Invariant Scattering for Texture Discrimination.  In CVPR, 2013.
  * [39] E. Simoncelli.  A rotation invariant pattern signature.  In ICIP, 1996.
  * [40] Sivic and Zisserman.  Video Google: a text retrieval approach to object matching in videos.  In ICCV, 2003.
  * [41] O. J. O. Söderkvist.  Computer Vision Classification of Leaves from Swedish Trees.  2001\.
  * [42] K. Sohn and H. Lee.  Learning Invariant Representations with Local Transformations.  In ICML, 2012.
  * [43] M. Weiler, F. A. Hamprecht, and M. Storath.  Learning Steerable Filters for Rotation Equivariant CNNs.  In CVPR, 2018.
  * [44] D. E. Worrall, S. J. Garbin, D. Turmukhambetov, and G. J. Brostow.  Harmonic Networks: Deep Translation and Rotation Equivariance.  In CVPR, 2017.