# Indirect predicates for geometric constructions

Marco Attene<EMAIL_ADDRESS>CNR-IMATI, Via De Marni 6, 16149
Genova, Italy

###### Abstract

Geometric predicates are a basic ingredient to implement a vast range of
algorithms in computational geometry. Modern implementations employ floating
point filtering techniques to combine efficiency and robustness, and state-of-
the-art predicates are guaranteed to be always exact while being only slightly
slower than corresponding (inexact) floating point implementations.
Unfortunately, if the input to these predicates is an intermediate
construction of an algorithm, its floating point representation may be
affected by an approximation error, and correctness is no longer guaranteed.

This paper introduces the concept of _indirect_ geometric predicate: instead
of taking the intermediate construction as an explicit input, an indirect
predicate considers the primitive geometric elements which are combined to
produce such a construction. This makes it possible to keep track of the
floating point approximation, and thus to exploit efficient filters and
expansion arithmetic to exactly resolve the predicate with minimal overhead
with respect to a naive floating point implementation. As a representative
example, we show how to extend standard predicates to the case of points of
intersection of linear elements (i.e. lines and planes) and show that, on
classical problems, this approach outperforms state-of-the-art solutions based
on lazy exact intermediate representations.

###### keywords:

Geometric predicates , Exact arithmetic , Filtering

##  1 Introduction

The study and design of computational geometry algorithms is widely based on
the so-called Real Arithmetic Model. According to that model, coordinates,
lengths, angles, and many other entities, are all considered real numbers, and
the algorithm correctness is based on this assumption. Practical
implementations, however, are typically based on floating point (FP)
arithmetic, and the aforementioned entities are represented by FP numbers.
Floating point operations produce approximated results and, in some cases, the
mismatch with the theoretical exact result might put the algorithm in an
inconsistent state, cause infinite loops, and lead to crashes [1] [2]. An easy
solution to this problem is to replace FP numbers with exact number types, but
the consequent slowdown is rarely acceptable for practical applications.
Geometric algorithms have the peculiar property that their flow depends on the
analysis of geometric entities and their configuration, and geometric
predicates are the basic building blocks to perform such an analysis. Broadly
speaking, a geometric predicate is a function that maps a set of geometric
entities to an element of a (small) set of values. For example, the predicate
collinear maps a triplet of points to a value in the set {true, false},
whereas the predicate orient2d maps a triplet of points to a value in the set
{left_turn, collinear, right_turn}. To avoid inconsistencies, it is often
essential that each predicate that determines the program flow produces exact
results.

Arithmetic filtering techniques (see Sect. 2.1) allow implementing predicates
so that they are always exact, while being only slightly slower than a
corresponding naive floating point implementation [3]. These techniques are
appropriate when the predicate parameters are all input values and hence not
computed as intermediate constructions by the algorithm. Conversely, if
intermediate constructions are used as predicate parameters, an approximation
during their computation may corrupt the result of the predicate,
independently of how cleverly it has been implemented. Unfortunately, many
useful algorithms require such intermediate constructions, and for these cases
state-of-the-art robust solutions mostly rely on exact or lazy exact
arithmetic kernels, which are still significantly slower and more memory
intensive than floating point implementations.

In this paper, we show that most of these algorithms can be robustly
implemented thanks to a novel concept of indirect geometric predicate. In
particular, we show that the geometric construction itself can be filtered
along with the predicate, thus leading to provably robust algorithms which are
comparable with naive floating point implementations in terms of both
performances and memory footprint. The reference implementation and the
scripts to reproduce the results in the paper are provided in the additional
material and will be released as an open-source project.

##  2 Background and prior art

The majority of geometric predicates evaluate the sign of a (homogeneous)
polynomial, typically a matrix determinant. As an example, the Cartesian 2D
orientation of three points $a$, $b$, and $c$ corresponds to the sign of the
determinant of the $2\times 2$ matrix whose rows are the vectors $b-a$ and
$c-a$. Instead of using exact arithmetic everywhere, which is too slow, modern
implementations use arithmetic filtering. These implementations evaluate the
polynomial sign using only FP arithmetic, but while keeping track of the
approximation error. If this error is large enough to make the sign ambiguous
(i.e. when the filter “fails”), they switch to exact arithmetic. The idea is
that the failure rate is low enough to make the impact of exact arithmetic
virtually negligible. Several approaches exist to bound the accumulated error.
In general, calculating a precise bound is slower but minimizes the
probability of filter failures. Conversely, conservatively approximating the
bound is faster but less precise.

###  2.1 Static, semi-static and dynamic filtering

Given an arithmetic expression involving a number of variables (e.g. the
polynomial of a predicate), its evaluation using FP arithmetic is subject to a
roundoff error that depends on the values assumed by the variables. This error
can be bounded by a constant threshold that, in turn, can be computed
analytically and initialized at compilation time. At runtime, if the magnitude
of the evaluated expression is larger than this threshold, its sign is
guaranteed to be correct. This approach is known as _static filtering_ [4],
and the overhead with respect to a plain FP implementation is due to a simple
comparison. Without any assumption on the input, the threshold is typically
too large and hence not really useful in practice. Nonetheless, if all the
potential input values are known at the beginning of the program, the
threshold can be computed once for all and then used to filter all the calls
to the predicate.

Since the actual input to a predicate is typically made of points which are
close to each other, considering only these points to derive a ”local”
threshold leads to tighter error bounds and, hence, to less filter failures.
This approach is known as _semi-static filtering_ , and the overhead is due to
a few additional operations to compute the threshold at each predicate call
[5].

Instead of computing an explicit threshold, interval arithmetic can be used to
tightly bound the exact result. Basically, the result of any single operation
is rounded by both excess and defect, thus leading to an interval containing
the exact result. All the operations accumulate these approximations into
growing intervals. If the interval representing the eventual result does not
contain the zero, then its sign is defined without ambiguity. This approach is
an example of _dynamic filtering_ , and evaluating the predicate in this way
takes from three to eight times longer than plain FP calculation [6].

It is often advantageous to use all these approaches in a cascaded multi-stage
evaluation. Thus, the first attempt is based on plain FP calculation, and a
single additional comparison allows implementing a static filter. If this
filter fails, a few more comparisons (e.g. four for the orient2d predicate)
allow implementing a semi-static filter. If this filter also fails, the
polynomial is re-evaluated using interval arithmetic, which is the last
attempt. As mentioned above, the first (static) filter makes sense only if a
reasonably tight bound can be computed based on the input. In general, the
process may safely start from the semi static filter while keeping the
overhead reasonably small.

###  2.2 Exact and lazy exact evaluation

When all the filters fail, the predicate expression must be evaluated using
exact arithmetic. State-of-the-art libraries allow, for example, to represent
intermediate results using rational numbers, thus without any loss of
precision [7]. Also, polynomials can be exactly evaluated by using arbitrary
precision floating point numbers [8] or so-called floating point expansions
[9]. This latter approach is particularly efficient as it fully exploits the
floating point hardware. Shewchuck’s geometric predicates [10] and Levy’s
Geogram [11] exploit floating point expansions and are nearly as efficient as
plain FP implementations on classical test beds. At the time of writing,
Geogram predicates can be considered to represent the state of the art.

When intermediate constructions are used as input to a predicate, current
robust solutions are based on lazy exact evaluation [12]. This is the approach
adopted in CGAL’s Lazy Exact kernel: any intermediate construction is
symbolically represented by the history of all the operations that combine
input values to build the construction [13]. Besides keeping track of the
operations, the construction is also approximately evaluated using interval
arithmetic [6]. Since predicate expressions are also evaluated using interval
arithmetic, their parameters can be intermediate constructions whose intervals
propagate throughout the predicate’s expression evaluation. If the resulting
interval contains the zero, it means that the sign is ambiguous, and all the
involved symbolic representations are tracked back and explicitly evaluated
using exact arithmetic. The expectation is that in most of the cases interval
arithmetic is sufficiently precise to reach an unambiguous conclusion, and
hence exact evaluation is necessary only in a minority of the cases.
Nonetheless, even in the ideal case where exact arithmetic is never called
into play, interval arithmetic can be up to eight times slower than FP
arithmetic [6], and keeping track of the whole history has its own cost too.

Note that all the mechanisms described so far have been designed to compensate
for numerical imprecision. Nonetheless, numerical imprecision is not the only
source of ambiguity in geometric computation: many algorithms assume that
their input is in so-called _general position_ and, if this is not the case,
this precondition needs to be simulated [14]. This is rather orthogonal to
compensating for numerical imprecision, and the two techniques should be
combined when necessary.

##  3 Geometric algorithm classification

Herewith, we classify geometric algorithms based on the number of steps that
require intermediate constructions. There are algorithms that have no such
steps, that is, where all the predicates depend directly on input values.
Then, there are algorithms having only one such step, where predicates depend
either on input values or on intermediate constructions that, in their turn,
depend directly on input values. And so on, up to arbitrarily complicated
algorithms with cascaded sequences of predicate evaluations and intermediate
constructions that mutually depend on each other. Our key observation is that
the vast majority of state-of-the-art practical algorithms belong to the first
two classes.

If $S$ is a set of 2D points whose coordinates are FP numbers, calculating
their Delaunay triangulation amounts to repeatedly evaluate orient2d and
incircle predicates whose input are points in $S$. Thus, this is a typical
example of algorithm of the first class. Other examples include calculating
convex hulls, Voronoi diagrams, surface interpolation, and many others.

If $S$ is a set of 2D segments that possibly intersect, calculating the
constrained Delaunay triangulation of their arrangement is an example of
algorithm of the second class. Indeed, the intersection of two segments is an
intermediate construction used as input by the orient2d and incircle
predicates. Other examples include the calculation of Minkowski sums [15],
mesh booleans [16], exact mesh repairing [17], and so on.

Current techniques allow combining robustness and efficiency only for
algorithms of the first class thanks to arithmetic filtering. In this paper,
we propose a novel approach to move forward and include algorithms of the
second class too, thus covering a much wider subset of fundamental geometric
algorithms.

##  4 Indirect geometric predicates

The core intuition behind our contribution is the following: if the code to
calculate the intermediate construction is included in the predicate itself,
arithmetic filtering may be applied to that code as well, and the predicate
can be evaluated using FP arithmetic as long as the filter does not fail.

When intermediate constructions are used as parameters for a predicate, they
are passed to the predicate in an unevaluated form: in this case, we say that
these parameters are _implicit_ and the predicate is _indirect_. An indirect
predicate must accept at least one implicit parameter, but may accept standard
_explicit_ parameters too. In the latter case, the parameter has an explicit
value and can be assumed to be exact (e.g. it is an input value or the result
of an operation that does not introduce approximation errors). Conversely, in
the former case the parameter is defined as an unevaluated function of
explicit values (e.g. the average of input values). Since most geometric
predicates have points as parameters, we extend this terminology to points and
say that an explicit point is represented by its explicit coordinates, whereas
an implicit point is represented by a geometric construction over explicit
points (e.g. by the intersection of two lines, each defined by two explicit
points).

The idea is to rewrite the expression of an implicit parameter $i$ as a
fraction $\lambda/d$ where both the numerator $\lambda$ and the denominator
$d$ are polynomials. Then, any standard predicate can be rewritten by
substituting, in its own polynomial $\Lambda$, explicit parameters with
polynomial fractions representing implicit parameters, thus reducing the whole
predicate expression to another polynomial fraction
$\Lambda^{\prime}/D^{\prime}$. Evaluating the sign of such a fraction amounts
to evaluate the sign of both the numerator $\Lambda^{\prime}$ and denominator
$D^{\prime}$, and both can be done robustly and efficiently through arithmetic
filtering.

When an implicit parameter is a point, we derive a common denominator for all
the fractions representing its Cartesian coordinates. Hence, an implicit point
$i=(i_{x},i_{y},i_{z})$ is rewritten as
$i=(\lambda_{ix}/d_{i},\lambda_{iy}/d_{i},\lambda_{iz}/d_{i})$ or, using
homogeneous coordinates, as $(\lambda_{ix},\lambda_{iy},\lambda_{iz},d_{i})$.

Note that possible singular/degenerate configurations must be detected and
appropriately handled. Typically, this amounts to check whether $d$ becomes
zero for each of the implicit input values.

In the following subsections, we show how to apply this transformation to some
non-trivial geometric constructions, and how to apply them to fundamental
geometric predicates.

###  4.1 Line-line intersections in 2D

Let p be an implicit 2D point defined as the intersection of two lines $l_{a}$
and $l_{b}$, each defined by two points. Specifically, $\textbf{e}_{a1}$ and
$\textbf{e}_{a2}$ define $l_{a}$, whereas $\textbf{e}_{b1}$ and
$\textbf{e}_{b2}$ define $l_{b}$ (Fig. 1-left). Therefore:

$\displaystyle\textbf{p}=(\frac{\lambda_{x}}{d},\frac{\lambda_{y}}{d})$

where:

$\displaystyle\lambda_{x}$
$\displaystyle=(e_{a1x}e_{a2y}-e_{a2x}e_{a1y})(e_{b1x}-e_{b2x})+$
$\displaystyle-(e_{b1x}e_{b2y}-e_{b2x}e_{b1y})(e_{a1x}-e_{a2x})$
$\displaystyle\lambda_{y}$
$\displaystyle=(e_{a1x}e_{a2y}-e_{a2x}e_{a1y})(e_{b1y}-e_{b2y})+$
$\displaystyle-(e_{b1x}e_{b2y}-e_{b2x}e_{b1y})(e_{a1y}-e_{a2y})$
$\displaystyle d$
$\displaystyle=(e_{a1x}-e_{a2x})(e_{b1y}-e_{b2y})-(e_{a1y}-e_{a2y})(e_{b1x}-e_{b2x})$

Note that p is undefined if $d=0$, which happens if $l_{a}$ and $l_{b}$ are
parallel or if either of the two point pairs is degenerate.

###  4.2 Line-plane intersections in 3D

Let p be an implicit 3D point defined as the intersection of a straight line
$L$ and a plane $\Pi$, let $L$ be defined by two points $\textbf{q}_{1}$,
$\textbf{q}_{2}$, and let $\Pi$ be defined by three points r, s, t in the
Euclidean space (Fig. 1-right). Therefore:

$\displaystyle\textbf{p}=(\frac{\lambda_{x}}{d},\frac{\lambda_{y}}{d},\frac{\lambda_{z}}{d})$

where

$\displaystyle d=\begin{vmatrix}\textbf{q}_{1}-\textbf{q}_{2}\\\
\textbf{s}-\textbf{r}\\\
\textbf{t}-\textbf{r}\end{vmatrix},n=\begin{vmatrix}\textbf{q}_{1}-\textbf{r}\\\
\textbf{s}-\textbf{r}\\\ \textbf{t}-\textbf{r}\end{vmatrix}$
$\displaystyle\lambda_{x}=dq_{1x}+nq_{2x}-nq_{1x}$
$\displaystyle\lambda_{y}=dq_{1y}+nq_{2y}-nq_{1y}$
$\displaystyle\lambda_{z}=dq_{1z}+nq_{2z}-nq_{1z}$

Note that p is undefined if $d=0$, which happens if $L$ and $\Pi$ are
parallel, if $\textbf{q}_{1}=\textbf{q}_{2}$, or if r, s, t are collinear or
coincident.

Figure 1: Implicit point examples. Left - line-line intersection: four
explicit 2D points define the implicit point p. Right - line-plane
intersection: five explicit 3D points define the implicit point p.

###  4.3 orient2d

Let $\textbf{p}_{1}$, $\textbf{p}_{2}$, $\textbf{p}_{3}$ be three points on
the Euclidean plane. The polynomial expression of the classical orient2d
predicate is
$\Lambda=det|\textbf{p}_{2}-\textbf{p}_{1};\textbf{p}_{3}-\textbf{p}_{1}|=(p_{2x}-p_{1x})(p_{3y}-p_{1y})-(p_{2y}-p_{1y})(p_{3x}-p_{1x})$.
To start with, consider the simple case where just one of the points is an
intermediate construction. Let it be $\textbf{p}_{1}$, and let us assume that
we are able to rewrite its Cartesian coordinates as polynomial fractions (e.g.
as described in Sect. 4.1). Therefore:

$\displaystyle p_{1x}=\frac{\lambda_{x}}{d},p_{1y}=\frac{\lambda_{y}}{d}$

By substituting $\textbf{p}_{1}$ with $(\lambda_{x}/d,\lambda_{y}/d)$, we
obtain
$\Lambda=(p_{2x}-\lambda_{x}/d)(p_{3y}-\lambda_{y}/d)-(p_{2y}-\lambda_{y}/d)(p_{3x}-\lambda_{x}/d)=\Lambda^{\prime}/D^{\prime}$,
where
$\Lambda^{\prime}=(dp_{2x}-\lambda_{x})(dp_{3y}-\lambda_{y})-(dp_{2y}-\lambda_{y})(dp_{3x}-\lambda_{x})$
and $D^{\prime}=d^{2}$. Since the denominator is a square, it can be safely
dropped and the sign is unambiguously defined by $\Lambda^{\prime}$.

Now, let us consider the most generic case where all the three points are
intermediate constructions. We obtain:

$\displaystyle
p_{1x}=\frac{\lambda_{1x}}{d_{1}},p_{1y}=\frac{\lambda_{1y}}{d_{1}}$
$\displaystyle
p_{2x}=\frac{\lambda_{2x}}{d_{2}},p_{2y}=\frac{\lambda_{2y}}{d_{2}}$
$\displaystyle
p_{3x}=\frac{\lambda_{3x}}{d_{3}},p_{3y}=\frac{\lambda_{3y}}{d_{3}}$

Where the $\lambda_{i}$ and the $d_{i}$ are computed as above. By substituting
$\textbf{p}_{1}$, $\textbf{p}_{2}$, $\textbf{p}_{3}$ with their corresponding
fractions, we obtain
$\Lambda=(\lambda_{2x}/d_{2}-\lambda_{1x}/d_{1})(\lambda_{3y}/d_{3}-\lambda_{1y}/d_{1})-(\lambda_{2y}/d_{2}-\lambda_{1y}/d_{1})(\lambda_{3x}/d_{3}-\lambda_{1x}/d_{1})=\Lambda^{\prime}/D^{\prime}$,
where

$\displaystyle\Lambda^{\prime}=(d_{1}\lambda_{2x}-d_{2}\lambda_{1x})(d_{1}\lambda_{3y}-d_{3}\lambda_{1y})-(d_{1}\lambda_{2y}-d_{2}\lambda_{1y})(d_{1}\lambda_{3x}-d_{3}\lambda_{1x})$
$\displaystyle D^{\prime}=d_{1}^{2}d_{2}d_{3}.$

###  4.4 incircle

Let $\textbf{p}_{1}$, $\textbf{p}_{2}$, $\textbf{p}_{3}$, $\textbf{p}_{4}$ be
four points on the Euclidean plane. The polynomial expression of the classical
incircle predicate is the following:

$\displaystyle\Lambda=\begin{vmatrix}p_{1x}-p_{4x}&p_{1y}-p_{4y}&(p_{1x}-p_{4x})^{2}+(p_{1y}-p_{4y})^{2}\\\
p_{2x}-p_{4x}&p_{2y}-p_{4y}&(p_{2x}-p_{4x})^{2}+(p_{2y}-p_{4y})^{2}\\\
p_{3x}-p_{4x}&p_{3y}-p_{4y}&(p_{3x}-p_{4x})^{2}+(p_{3y}-p_{4y})^{2}\end{vmatrix}$

If each of the four points is implicitly defined, we can rewrite this
expression as:

$\displaystyle\Lambda=\begin{vmatrix}\frac{\lambda_{1x}}{d_{1}}-\frac{\lambda_{4x}}{d_{4}}&\frac{\lambda_{1y}}{d_{1}}-\frac{\lambda_{4y}}{d_{4}}&(\frac{\lambda_{1x}}{d_{1}}-\frac{\lambda_{4x}}{d_{4}})^{2}+(\frac{\lambda_{1y}}{d_{1}}-\frac{\lambda_{4y}}{d_{4}})^{2}\\\
\frac{\lambda_{2x}}{d_{2}}-\frac{\lambda_{4x}}{d_{4}}&\frac{\lambda_{2y}}{d_{2}}-\frac{\lambda_{4y}}{d_{4}}&(\frac{\lambda_{2x}}{d_{2}}-\frac{\lambda_{4x}}{d_{4}})^{2}+(\frac{\lambda_{2y}}{d_{2}}-\frac{\lambda_{4y}}{d_{4}})^{2}\\\
\frac{\lambda_{3x}}{d_{3}}-\frac{\lambda_{4x}}{d_{4}}&\frac{\lambda_{3y}}{d_{3}}-\frac{\lambda_{4y}}{d_{4}}&(\frac{\lambda_{3x}}{d_{3}}-\frac{\lambda_{4x}}{d_{4}})^{2}+(\frac{\lambda_{3y}}{d_{3}}-\frac{\lambda_{4y}}{d_{4}})^{2}\end{vmatrix}$

The common denominator to all the nine elements is
$D^{\prime}=(d_{1}d_{2}d_{3}d_{4})^{2}$. When none of the $d$’s is zero, the
sign of $\Lambda$ corresponds to the sign of the following $\Lambda^{\prime}$:

$\displaystyle\Lambda^{\prime}=\begin{vmatrix}m_{1,1}&m_{1,2}&m_{1,3}\\\
m_{2,1}&m_{2,2}&m_{2,3}\\\ m_{3,1}&m_{3,2}&m_{3,3}\end{vmatrix}$
$\displaystyle m_{1,1}$
$\displaystyle=d_{1}d_{4}^{2}\lambda_{1x}-d_{1}^{2}d_{4}\lambda_{4x}$
$\displaystyle m_{1,2}$
$\displaystyle=d_{1}d_{4}^{2}\lambda_{1y}-d_{1}^{2}d_{4}\lambda_{4y}$
$\displaystyle m_{1,3}$
$\displaystyle=d_{4}^{2}(\lambda_{1x}^{2}+\lambda_{1y}^{2})+d_{1}^{2}(\lambda_{4x}^{2}+\lambda_{4y}^{2})-2d_{1}d_{4}(\lambda_{1x}\lambda_{4x}+\lambda_{1y}\lambda_{4y})$
$\displaystyle m_{2,1}$
$\displaystyle=d_{2}d_{4}^{2}\lambda_{2x}-d_{2}^{2}d_{4}\lambda_{4x}$
$\displaystyle m_{2,2}$
$\displaystyle=d_{2}d_{4}^{2}\lambda_{2y}-d_{2}^{2}d_{4}\lambda_{4y}$
$\displaystyle m_{2,3}$
$\displaystyle=d_{4}^{2}(\lambda_{2x}^{2}+\lambda_{2y}^{2})+d_{2}^{2}(\lambda_{4x}^{2}+\lambda_{4y}^{2})-2d_{2}d_{4}(\lambda_{2x}\lambda_{4x}+\lambda_{2y}\lambda_{4y})$
$\displaystyle m_{3,1}$
$\displaystyle=d_{3}d_{4}^{2}\lambda_{3x}-d_{3}^{2}d_{4}\lambda_{4x}$
$\displaystyle m_{3,2}$
$\displaystyle=d_{3}d_{4}^{2}\lambda_{3y}-d_{3}^{2}d_{4}\lambda_{4y}$
$\displaystyle m_{3,3}$
$\displaystyle=d_{4}^{2}(\lambda_{3x}^{2}+\lambda_{3y}^{2})+d_{3}^{2}(\lambda_{4x}^{2}+\lambda_{4y}^{2})-2d_{3}d_{4}(\lambda_{3x}\lambda_{4x}+\lambda_{3y}\lambda_{4y})$

###  4.5 orient2d3d

Let $\Pi$ be an oriented plane identified by three non-aligned points r, s and
t in the Euclidean 3D space. Also, let $\textbf{p}_{1}$, $\textbf{p}_{2}$,
$\textbf{p}_{3}$ be three points on $\Pi$. The orient2d3d predicate
generalizes the standard orient2d predicate, and states whether
$\textbf{p}_{1}$, $\textbf{p}_{2}$, $\textbf{p}_{3}$ form a right turn, a left
turn, or a straight line on $\Pi$. It turns out that
$\texttt{orient2d3d}(\textbf{p}_{1},\textbf{p}_{2},\textbf{p}_{3},\textbf{r},\textbf{s},\textbf{t})=\texttt{orient2d}(\Gamma(\textbf{p}_{1}),\Gamma(\textbf{p}_{2}),\Gamma(\textbf{p}_{3}))*\texttt{orient2d}(\Gamma(\textbf{r}),\Gamma(\textbf{s}),\Gamma(\textbf{t}))$,
where $\Gamma(\textbf{p}_{i})$ indicates the 2D projection of $\textbf{p}_{i}$
on one of the three axis-aligned planes such that
$\texttt{orient2d}(\Gamma(\textbf{r}),\Gamma(\textbf{s}),\Gamma(\textbf{t}))\neq
0$. An effective way to select an appropriate projection plane is to consider
a normal vector at $\Pi$ and simply drop the coordinate with largest
magnitude. Without loss of generality, from now on we assume that we are
dropping the Z coordinate or, equivalently, that we are projecting on the XY
plane (i.e. $\Gamma(<p_{x},p_{y},p_{z}>)=<p_{x},p_{y}>$).

Now, let us suppose that $\textbf{p}_{1}$ is implicitly defined (e.g. as
described in Sect. 4.2). Therefore:

$\displaystyle
p_{1x}=\frac{\lambda_{x}}{d},p_{1y}=\frac{\lambda_{y}}{d},p_{1z}=\frac{\lambda_{z}}{d}$

Since we are projecting on XY, we do not even need to calculate $\lambda_{z}$,
and simply re-write
$\Lambda=(p_{2x}-p_{1x})(p_{3y}-p_{1y})-(p_{2y}-p_{1y})(p_{3x}-p_{1x})=(p_{2x}-\lambda_{x}/d)(p_{3y}-\lambda_{y}/d)-(p_{2y}-\lambda_{y}/d)(p_{3x}-\lambda_{x}/d)=\Lambda^{\prime}/D^{\prime}$,
where
$\Lambda^{\prime}=(dp_{2x}-\lambda_{x})(dp_{3y}-\lambda_{y})-(dp_{2y}-\lambda_{y})(dp_{3x}-\lambda_{x})$
and $D^{\prime}=d^{2}$. Hence,
$\texttt{orient2d}(\Gamma(\textbf{p}_{1}),\Gamma(\textbf{p}_{2}),\Gamma(\textbf{p}_{3}))=\texttt{sign}(\Lambda^{\prime})$.

###  4.6 orient3d

Similar arguments hold for the orient3d predicate. When all its input points
are explicit, the orient3d polynomial is:

$\displaystyle\Lambda=\begin{vmatrix}p_{1x}-p_{4x}&p_{1y}-p_{4y}&p_{1z}-p_{4z}\\\
p_{2x}-p_{4x}&p_{2y}-p_{4y}&p_{2z}-p_{4z}\\\
p_{3x}-p_{4x}&p_{3y}-p_{4y}&p_{3z}-p_{4z}\end{vmatrix}$

If we replace $\textbf{p}_{1}$ with an implicit point we obtain the following
expression for an indirect orient3d predicate:

$\displaystyle\Lambda^{\prime\prime}=\begin{vmatrix}\lambda_{x}-dp_{4x}&\lambda_{y}-dp_{4y}&\lambda_{z}-dp_{4z}\\\
dp_{2x}-dp_{4x}&dp_{2y}-dp_{4y}&dp_{2z}-dp_{4z}\\\
dp_{3x}-dp_{4x}&dp_{3y}-dp_{4y}&dp_{3z}-dp_{4z}\end{vmatrix}$

which, if $d\neq 0$, has the same sign as the simpler:

$\displaystyle\Lambda^{\prime}=\begin{vmatrix}\lambda_{x}-dp_{4x}&\lambda_{y}-dp_{4y}&\lambda_{z}-dp_{4z}\\\
p_{2x}-p_{4x}&p_{2y}-p_{4y}&p_{2z}-p_{4z}\\\
p_{3x}-p_{4x}&p_{3y}-p_{4y}&p_{3z}-p_{4z}\end{vmatrix}$

where the corresponding $D^{\prime}$ is $d$.

Formulas for the case where more than one point is implicit can be easily
derived using a similar rewriting.

##  5 Computational models and implementation

The input to an indirect predicate is a set of points, each being either
explicit or implicit. In the former case, the point is a standard array of
Cartesian coordinates, each encoded as a double precision FP number; in the
latter case, it is a collection of explicit points to be combined. Upon
invocation, an indirect predicate determines which of its input parameters are
implicit, and internally calculates $\lambda$ and $d$ for each of them. In our
implementation we consider three computation models: floating point arithmetic
(with semi-static filter), interval arithmetic (with dynamic filter),
expansion arithmetic (exact). The predicate is evaluated with the fastest
model which guarantees exactness.

###  5.1 Filtered evaluation in floating point

We begin by calculating the $d$’s and, for each of them, use a semi-static
filter that states whether its sign is unambiguous. If the filter succeeds for
all the $d$’s, we compute the $\lambda$’s, evaluate $\Lambda^{\prime}$, and
use another semi-static filter to guarantee that its sign is correct. If the
former filter fails for one of the $d$’s, it means that precision is
insufficient to exclude that the implicit point is undefined (e.g.
intersection of parallel lines). In this case, it makes no sense to proceed
with this arithmetic model. As soon as a filter fails, we stop and switch to
interval arithmetic (Sect. 5.2). Upon success, the sign of $\Lambda^{\prime}$
is possibly inverted depending on the sign of $D^{\prime}$. We observe that
$D^{\prime}$ is just a product of the $d$’s, therefore we simply count how
many of them have a negative sign (possible multiplicities must be
considered). If they are an even number, we return the sign of
$\Lambda^{\prime}$, otherwise its inverse. Semi-static filters for the
floating point model are calculated using a variation of [5] (see Appendix A
for details).

###  5.2 Filtered evaluation using intervals

When starting with this arithmetic model, we first set the CPU rounding mode
to towards_infinity. This is required to correctly and efficiently evaluate
the subsequent interval arithmetic block. Then, we calculate the $d$’s and
check that their sign is not ambiguous (i.e. their intervals do not contain
the zero). If this dynamic filter succeeds for all the $d$’s, we compute the
$\lambda$’s, evaluate $\Lambda^{\prime}$, and check that the resulting sign is
not ambiguous. As soon as one of the filters fail we stop and switch to exact
arithmetic (Sect. 5.3). If all the filters succeed, the sign of
$\Lambda^{\prime}$ is possibly inverted depending on the sign of the $d$’s.
Upon termination (both with success or failure), we reset the rounding to the
standard to_nearest mode required by the other computational models.

###  5.3 Exact evaluation using expansions

When the predicate is computed using an exact arithmetic model there is a
minor difference. Indeed, when computing the $d$’s it might happen that some
of them are exactly zero. If so, the predicate cannot be evaluated because at
least one of its implicit parameters is undefined. Conversely, if all the
$d$’s are nonzero, we calculate $\Lambda^{\prime}$ and check whether it is
zero. If so, there is no need to check the $d$’s and the predicate returns
zero. Otherwise, we possibly invert the sign as in the previous models.

Since $\Lambda^{\prime}$ is a polynomial, we may exactly evaluate its sign
using floating point expansions. Roughly speaking, an expansion is an array of
FP numbers whose sum is the exact number being represented. Under certain
assumptions, FP expansions can be efficiently summed, subtracted, and
multiplied, all without error and while exploiting the floating point hardware
[9].

###  5.4 Caching

If the same implicit point is used as input to many predicate calls, its
$\lambda$’s and $d$ are recomputed many times. Therefore, upon the first
occurrence it is worth storing their values for future reuse. For the sake of
implementing filters, besides these values we must also store information to
accumulate error bounds. For semi-static filters, this amounts to store an
additional FP number (see Appendix A). For dynamic filters, this amounts to
store the $\lambda$’s and $d$ using intervals.

Clearly, such a caching mechanism speeds up the computation at the expense of
a higher memory usage. Our experiments (see Sect. 7) show that caching
$\lambda$’s and $d$’s for floating point and interval arithmetic models
produces relevant speed ups while requiring a relatively low increase in
memory footprint. Conversely, storing these values for the exact model is not
advantageous: indeed, storing expansions for each point requires a significant
amount of memory, and the time overhead due to the management of all this
additional data is higher than what we save in the few switches which are
required in typical cases (Table 3). Note that pre-calculated intervals are
also used in CGAL Lazy exact kernel [18], but the use of cached FP filters
makes our predicate evaluaton significantly faster (Sect. 7).

##  6 Code generation

An indirect predicate must consider all the possible configurations of its
input parameters. For example, the orient2d predicate has three points as
parameters, and the type of each of these points can be either _implicit_ or
_explicit_. The most generic implementation may assume that all the points are
implicit and, if any of them is explicit, its $\lambda$’s can be simply
replaced with the Cartesian coordinates and its $d$ with 1. Though this
approach is correct, such a naive implementation would lead to unnecessarily
loose filters which, in turn, would require unnecessary switches to more
precise models. Consequently, this implementation would be unnecessarily slow.

A much more efficient approach is based on the implementation of a specific
instance of the predicate for each of the possible input configurations. If
the three input points of our indirect orient2d predicate are all explicit,
the predicate reduces to the standard orient2d, which is very efficient and
accumulates small errors. If only the first point is implicit, we may use a
dedicated implementation based on the $\Lambda$ and $D$ described earlier in
Sect. 4.3. This is slightly less efficient than the previous version and
produces a slightly higher error. And so on. Hence, an indirect orient2d
predicate might be efficiently implemented using eight instances (EEE, EEI,
EIE, EII, IEE, IEI, IIE, III, where the three characters indicate the type of
each input point). Actually, many of these instances can be reduced to each
other by transposing their input parameters and possibly inverting the
resulting sign (e.g. orient2d(a, b, c) = orient2d(b, c, a) = -orient2d(a, c,
b)). It is easy to verify that four versions (EEE, IEE, IIE, III) of orient2d
are sufficient to cover all the cases.

Implementing even a single instance with proper semi-static and dynamic
filtering, and with expansion-based evaluation, is a tedious and error-prone
process. Fortunately, this process can be mostly automated [19]. If the
predicate is naively implemented as a C function using FP arithmetic, software
tools exist [5] to parse such a function and produce another C function which
includes a semi-static filter. These tools analyze the polynomial and
calculate an upper bound for the roundoff error. In the produced code, the
magnitude of the resulting determinant is compared with such a bound and, if
smaller, the predicate returns a particular value indicating that the filter
failed (e.g. FPG_UNCERTAIN). Meyer and Pion’s tool has been used and extended
by Levy in the _Predicate Construction Kit_ [11], where both the semi-static
filtered and the exact versions of the predicate are produced.

We have further extended upon this approach, and have implemented a software
that parses a text-based formula representing the polynomial. Starting from
such a formula, our tool produces three C++ functions, each implementing the
formula using one of the arithmetic models described earlier. A further fourth
function is generated that calls the other three in order and terminates as
soon as one of them returns a guaranteed result.

##  7 Experiments

For the sake of experimentation, we have implemented a C++ library with an API
defining a genericPoint class. A genericPoint can be specialized to be either
2D or 3D, implicit or explicit. The library also provides functions
implementing both standard and indirect predicates. Each such function takes a
number of genericPoint objects as input, internally determines their type, and
properly handles all the cases as described in Sect. 6. The library has been
used as a basis to implement a simple Delaunay triangulation algorithm that
can deal with genericPoint objects. Both the library and the triangulation
algorithm have been compiled and tested on an MS Windows 10 OS using Visual
C++ 2019. The objectives are measuring the efficiency of our predicates, the
overhead of indirect predicates when compared to their direct versions, the
impact of all the possible caching levels, and the improvement with respect to
the state of the art. Tests were run on a standard PC equipped with an Intel
Core i7-4770 with 16 Gb RAM. Based on this setting, the following experiments
were run. Results are summarized in Tables 1 and 2.

N. Points | Exp. 1.1 | Exp. 2.1
---|---|---
1000 | 0.000691 (2.74) | 0.000862 (2.77)
10000 | 0.008698 (6.29) | 0.008864 (6.49)
100000 | 0.117321 (27.99) | 0.12001 (25.56)
1000000 | 1.79624 (268.53) | 1.82953 (284.55)
Table 1: Results of our experiments 1.1 and 2.1. Elapsed time in seconds and, in parenthesis, peak memory usage in Mbytes. % implicit | Exp. 1.2 | Exp. 1.3 | Exp. 2.2 | Exp. 2.3
---|---|---|---|---
0 | 1.79624 (268.53) | 2.63308 (268.36) | 1.82953 (284.55) | 2.64815 (284.37)
10 | 2.23134 (292.84) | 5.0791 (292.99) | 2.56663 (321.96) | 7.14617 (322.95)
25 | 2.82864 (328.84) | 8.45898 (329.3) | 3.41211 (377.43) | 16.6439 (379.75)
50 | 3.62846 (389.14) | 14.0788 (390.61) | 4.24588 (469.68) | 38.3586 (472.40)
100 | 3.81999 (509.93) | 27.1925 (511.29) | 4.29742 (654.52) | 75.0639 (657.06)
Table 2: Results of our experiments 1.2, 1.3, 2.2 and 2.3. Elapsed time in
seconds and, in parenthesis, peak memory usage in Mbytes.

###  7.1 Experiment 1.1 - Random explicit 2D points

For this experiment, we created and triangulated a set of genericPoint
objects, each specialized to represent an explicit 2D point. Their Cartesian
coordinates are randomly sampled within the unit square $[0,1]\times[0,1]$.
Since all the points are explicit, the orient2d and incircle predicates used
by the algorithm always switch to the standard _direct_ version of the
predicates. A parameter allows controlling the cardinality of the point set in
the range [1, 1000000]. The objective of this test is to define a reference
for the other tests, in order to evaluate the overhead due to the use of
indirect predicates when compared with their classical direct version.

###  7.2 Experiment 1.2 - Random mixed 2D points

Herewith, our set of genericPoint objects still contains 2D points, but their
type is mixed. A subset of them are explicit, whereas the others are
implicitly defined as the intersection of line pairs. Explicit points are
randomly scattered within the unit square as in the previous experiment. Also,
line pairs are created so that their intersections are randomly scattered
within the unit square. A parameter allows controlling the percentage of
implicit points over a total of 1 million points constituting the entire set.
In this case, both orient2d and incircle may either switch to the direct
version or select the proper indirect predicate instance, depending on the
type of the genericPoints.

###  7.3 Experiment 1.3 - Regular mixed 2D points

This dataset of 1000000 genericPoint objects contains mixed 2D points as in
Experiment 1.2 but, instead of being randomly scattered, they form a regular
grid within the unit square. A subset of them are explicit, whereas the others
are implicitly defined as the intersection of line pairs. Even in this case, a
parameter allows controlling the percentage of implicit points over the total.
This dataset stresses the predicates that are largely required to switch to
exact arithmetic.

###  7.4 Experiment 2.1 - Random explicit 3D points

For this experiment, we created a set of genericPoint objects, each
specialized to represent an explicit 3D point whose Cartesian coordinates are
randomly sampled within the unit cube $[0,1]\times[0,1]\times[0,1]$. The
Delaunay triangulation is calculated on the projection of these points on the
XY plane. As for Experiment 1.1, the objective is to define a reference for
the following Experiments 2.2 and 2.3.

###  7.5 Experiment 2.2 - Random mixed 3D points

Here, our set of genericPoint objects contains 3D points, both explicit and
implicit, and the Delaunay triangulation is calculated on their projection on
the XY plane. As in the previous experiment, the first half of the set is made
of explicit points, whereas the other half are implicitly defined as the
intersection of lines and planes. Explicit points are randomly scattered
within the unit cube $[0,1]\times[0,1]\times[0,1]$. Line-plane pairs are
created so that their intersections belong to the unit cube. A parameter
allows controlling the percentage of implicit points over a total of 1 million
points constituting the entire set.

###  7.6 Experiment 2.3 - Regular mixed 3D points

This dataset of 1000000 genericPoint objects contains mixed 3D points as in
Experiment 2.2 but, instead of being randomly scattered, their projections on
the XY plane form a regular grid within the unit square. Their Z coordinate is
still randomnly selected in the interval [0,1]. A subset of the points are
explicit, whereas the others are implicitly defined as the intersection of
properly configured line-plane pairs. Even in this case, a parameter allows
controlling the percentage of implicit points over the total. This dataset
stresses the predicates that are largely required to switch to exact
arithmetic.

% implicit | no cache | FP | Interval | Exact
---|---|---|---|---
0 | 1.79552 (268.51) | 1.78396 (268.53) | 1.79624 (268.53) | 1.82953 (268.53)
10 | 2.49302 (283.07) | 2.4852 (286.20) | 2.23134 (292.84) | 2.27749 (358.66)
25 | 3.5918 (304.77) | 3.39718 (312.93) | 2.82864 (328.84) | 2.86765 (494.35)
50 | 5.25428 (341.23) | 5.00143 (357.23) | 3.62846 (389.14) | 3.79824 614.96)
100 | 7.16629 (413.65) | 6.9816 (445.84) | 3.81999 (509.93) | 3.9589 (1167.84)
Table 3: Results of our experiment 1.2 using different caching approaches.
$\lambda$’s and $d$ were cached only for the FP model (column FP), for both
the FP and interval arithmetic models (column Interval), or for all the
available models (column Exact). Elapsed time in seconds and, in parenthesis,
peak memory usage in Mbytes.

###  7.7 Comparison with state of the art

CGAL [18] is probably one of the most used libraries when robustness and
efficiency must be combined, in particular when important fundamental
algorithms come into play (e.g. Delaunay triangulations). That is why we
compare our results with CGAL. In our comparison, we always strive to exploit
CGAL in its most efficient configuration while providing a guarantee of
correctness. Namely, we want to guarantee that the resulting Delaunay
triangulation has the correct topology.

We re-run all our tests using CGAL Delaunay triangulation algorithm. For
experiments 1.1 and 2.1, where all the points are explicit, we configure CGAL
to use the Exact_Predicates_Inexact_Constructions kernel, whereas for the
others we need exact constructions to represent the intersections. Implicit
points are generated using CGAL::intersection(). For experiments that require
projection on the XY plane, CGAL::Projection_traits_xy_3 was used. We assume
that CGAL implements efficient state-of-the-art algorithms but, since using
two different implementations may introduce a bias independently of our
contribution, we have also re-implemented our own triangulation algorithm
using CGAL’s exact number type CGAL::Lazy_exact_nt<CGAL::Gmpq> and its
predicates CGAL::orientation() and CGAL::side_of_oriented_circle(). However,
after having verified that this re-implementation is slower than CGAL’s native
algorithm, we have abandoned it. Tables 4 and 5 report the results of our
experiments using CGAL.

N. Points | Exp. 1.1 | Exp. 2.1
---|---|---
1000 | 0.000253 (3.52) | 0.000273 (3.29)
10000 | 0.015636 (4.68) | 0.015622 (4.72)
100000 | 0.098743 (16.64) | 0.109367 (19.62)
1000000 | 1.168691 (135.81) | 1.15617 (167.65)
Table 4: Results of our experiments 1.1 and 2.1 using CGAL with exact predicates and inexact constructions kernel. Elapsed time in seconds and, in parenthesis, peak memory usage in Mbytes. % implicit | Exp. 1.2 | Exp. 1.3 | Exp. 2.2 | Exp. 2.3
---|---|---|---|---
0 | 6.92141 (168.86) | 26.4514 (362.26) | 33.2009 (193.40) | 68.7454 (483.43)
10 | 8.09864 (252.48) | 27.9512 (492.32) | 35.0445 (338.69) | 76.5043 (716.04)
25 | 9.76497 (377.75) | 32.0291 (698.53) | 38.7077 (553.73) | 88.9425 (1060.50)
50 | 11.9211 (586.53) | 40.966 (1040.29) | 42.2784 (908.59) | 111.112 (1634.54)
100 | 15.9442 (1006.73) | 60.9706 (1682.76) | 50.1841 (1605.84) | crashed
Table 5: Results of our experiments 1.2, 1.3, 2.2 and 2.3 using CGAL with
exact predicates and exact constructions kernel. Elapsed time in seconds and,
in parenthesis, peak memory usage in Mbytes.

###  7.8 Results discussion

In all the experiments, the total number of predicate calls depends on the
input size. For Experiment 1.1, these numbers are as follows: 1000 points
$\Rightarrow$ 4395 calls to orient2d and 7502 calls to incircle; 10K points
$\Rightarrow$ 45K and 94K calls; 100K points $\Rightarrow$ 469K and 1.14M
calls; 1M points $\Rightarrow$ 4.8M and 13.3M calls. These numbers do not vary
significantly across different experiments with same input size.

Tables 1 and 2 show that the impact of having implicit points in the input is
rather limited. The maximum increase in both time and memory usage is around
2x in 2D and 2.5x in 3D, which happens when all the points are implicit. Most
importantly, the resource demand gradually increases as the relative amount of
implicit points in the set grows. Hence, our approach is particularly
convenient in cases where implicit points are sparse. Experiments 1.3 and 2.3
show that our predicates behave reasonably well even in extremely degenerate
cases where exact arithmetic is used extensively. In particular, we observe
that the peak memory usage in these cases is virtually equivalent to that of
corresponding non-degenerate input sets. The impact on the elapsed time is of
course more significant but, even in this case, we observe that it gracefully
grows as the percentage of implicit points increases. Table 3 shows the impact
of caching the $\lambda$’s and $d$’s: when these values are stored for the FP
model only the overall gain in speed is modest, but when the Interval model is
included the algorithm becomes nearly twice as fast with only 23% more memory
requirements.

Our results also demonstrate that indirect predicates outperform state of the
art methods with same guarantees. Even if CGAL is faster when inexact
constructions are allowed (Experiments 1.1 and 2.1, see Tables 1 and 4), its
lazy exact kernel becomes necessary to represent implicit points and, in this
case, our implementation is from 3.86 to 11.7 times faster than CGAL, while
requiring up to 3.5 times less memory (see Tables 2 and 5 and Fig. 2).

Figure 2: Elapsed time for experiments 1.2, 1.3, 2.2 and 2.3 using indirect
predicates (blue) and CGAL (orange).

##  8 Conclusions

We have shown that, even if intermediate geometric constructions are
necessary, the efficiency of floating point hardware can be exploited to
quickly evaluate geometric predicates with guarantees of correctness. This
partly solves a problem introduced by Meyer and Pion [5] more than ten years
ago. The problem is only _partly_ solved because our approach requires that
intermediate constructions can be represented as polynomial fractions. Though
this is sufficient for most of the algorithms of the second type introduced in
Sect. 3, there are still a few useful algorithms that do not satisfy this
requirement. One noticeable example is the constrained Delaunay triangulation
in 3D, where state of the art algorithms use irrational numbers to generate
Steiner points [20]. Furthermore, more _numerically complex_ techniques exist
which involve cascaded predicates and constructions, and hence do not belong
to our second class of algorithms (e.g. [21]). Unfortunately, indirect
predicates do not help in these cases, and attempting to generalize their
paradigm is a challenging research direction: indeed, the combination of
cascaded constructions would easily lead to rather complicated expressions
whose filters become too large to be useful in practice. Nonetheless, indirect
predicates enable much faster implementations for a wide range of useful
algorithms of the second class.

We also observe that indirect predicates allow implementing algorithms with
guarantees of robustness and combinatorial correctness. Nonetheless, if their
results are the input to other algorithms which expect floating point
coordinates, implicit points must be evaluated and their coordinates
approximated. Such an approximation may lead to invalid configurations (e.g.
with self intersections) even if the combinatorics is exact. This problem is
known as _snap rounding_ and, though effective algorithms exist for the 2D
case, no practical guaranteed implementation exists yet for the 3D case. This
problem is worth a further investigation [22].

## Acknowledgements

This work is partly supported by the EU ERC Advanced Grant CHANGE, agreement
No 694515. Thanks are due to the colleagues at CNR-IMATI for helpful
discussions.

## References

## References

  * [1] L. Kettner, K. Mehlhorn, S. Pion, S. Schirra, C. Yap, Classroom examples of robustness problems in geometric computations, in: S. Albers, T. Radzik (Eds.), Algorithms – ESA 2004, Springer Berlin Heidelberg, Berlin, Heidelberg, 2004, pp. 702–713.
  * [2] C. Li, S. Pion, C. Yap, Recent progress in exact geometric computation, The Journal of Logic and Algebraic Programming 64 (1) (2005) 85 – 111, practical development of exact real number computation.  doi:https://doi.org/10.1016/j.jlap.2004.07.006.
  * [3] O. Devillers, F. P. Preparata, A probabilistic analysis of the power of arithmetic filters, Discrete & Computational Geometry 20 (4) (1998) 523–547.  doi:10.1007/PL00009400.
  * [4] S. Fortune, C. J. Van Wyk, Static analysis yields efficient exact integer arithmetic for computational geometry, ACM Trans. Graph. 15 (3) (1996) 223–248.  doi:10.1145/231731.231735.
  * [5] A. Meyer, S. Pion, FPG: A code generator for fast and certified geometric predicates, in: Real Numbers and Computers, 2008, pp. 47–60.
  * [6] H. Brönnimann, C. Burnikel, S. Pion, Interval arithmetic yields efficient dynamic filters for computational geometry, in: Proceedings of the Fourteenth Annual Symposium on Computational Geometry, SCG ’98, ACM, New York, NY, USA, 1998, pp. 165–174.  doi:10.1145/276884.276903.
  * [7] S. Fortune, C. J. Van Wyk, Efficient exact arithmetic for computational geometry, in: Proceedings of the Ninth Annual Symposium on Computational Geometry, SCG ’93, ACM, New York, NY, USA, 1993, pp. 163–172.  doi:10.1145/160985.161015.
  * [8] L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, P. Zimmermann, Mpfr: A multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Softw. 33 (2).  doi:10.1145/1236463.1236468.
  * [9] M. Joldes, O. Marty, J.-M. Muller, V. Popescu, Arithmetic algorithms for extended precision using floating-point expansions, IEEE TRANSACTIONS ON COMPUTERS 65 (4) (2016) 1197–1210.
  * [10] J. R. Shewchuk, Adaptive precision floating-point arithmetic and fast robust geometric predicates, Discrete & Computational Geometry 18 (3) (1997) 305–363.
  * [11] B. Lévy, Robustness and efficiency of geometric programs, Comput. Aided Des. 72 (C) (2016) 3–12.  doi:10.1016/j.cad.2015.10.004.
  * [12] S. Pion, A. Fabri, A generic lazy evaluation scheme for exact geometric computations, Science of Computer Programming 76 (4) (2011) 307 – 323, special issue on library-centric software design (LCSD 2006).  doi:https://doi.org/10.1016/j.scico.2010.09.003.
  * [13] S. Schirra, On the use of adaptive, exact decisions number types based on expression-dags in geometric computing, in: Proceedings of the 26th Canadian Conference on Computational Geometry, CCCG ’14, 2014.
  * [14] H. Edelsbrunner, E. P. Mücke, Simulation of simplicity: A technique to cope with degenerate cases in geometric algorithms, ACM Trans. Graph. 9 (1) (1990) 66–104.  doi:10.1145/77635.77639.
URL https://doi.org/10.1145/77635.77639

  * [15] P. Hachenberger, Exact minkowksi sums of polyhedra and exact and efficient decomposition of polyhedra into convex pieces, Algorithmica 55 (2) (2009) 329–345.  doi:10.1007/s00453-008-9219-6.
  * [16] Q. Zhou, E. Grinspun, D. Zorin, A. Jacobson, Mesh arrangements for solid geometry, ACM Trans. Graph. 35 (4) (2016) 39:1–39:15.  doi:10.1145/2897824.2925901.
  * [17] M. Attene, As-exact-as-possible repair of unprintable stl files, Rapid Prototyping Journal 24 (5) (2018) 855–864.  doi:10.1108/RPJ-11-2016-0185.
  * [18] The CGAL Project, CGAL User and Reference Manual, 5.0 Edition, CGAL Editorial Board, 2019.
URL https://doc.cgal.org/5.0/Manual/packages.html

  * [19] A. Nanevski, G. Blelloch, R. Harper, Automatic generation of staged geometric predicates, Higher-Order and Symbolic Computation 16 (4) (2003) 379–400.  doi:10.1023/A:1025876920522.
  * [20] H. Si, Tetgen, a delaunay-based quality tetrahedral mesh generator, ACM Trans. Math. Softw. 41 (2) (2015) 11:1–11:36.  doi:10.1145/2629697.
  * [21] F. Cazals, A. Parameswaran, S. Pion, Robust construction of the three-dimensional flow complex, in: Proceedings of the Twenty-Fourth Annual Symposium on Computational Geometry, SCG ’08, Association for Computing Machinery, New York, NY, USA, 2008, p. 182–191.  doi:10.1145/1377676.1377705.
URL https://doi.org/10.1145/1377676.1377705

  * [22] O. Devillers, S. Lazard, W. J. Lenhart, 3D Snap Rounding, in: B. Speckmann, C. D. Tóth (Eds.), 34th International Symposium on Computational Geometry (SoCG 2018), Vol. 99 of Leibniz International Proceedings in Informatics (LIPIcs), Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany, 2018, pp. 30:1–30:14.  doi:10.4230/LIPIcs.SoCG.2018.30.

##  Appendix A Semi-static filter calculation

A semi-static filter needs to (1) calculate a threshold value as a function of
input parameters $v_{i}$ and (2) verify that $|\Lambda^{\prime}|$ is larger
than this value. In our code, the threshold $\epsilon(v_{1},...,v_{n})$ is
calculated as a product $\delta(1)B(v_{1},...,v_{n})$, where $\delta(1)$ is a
constant value which depends only on the predicate’s expression and is
computed at compile time as described in [5]. Conversely, $B(v_{1},...,v_{n})$
depends on the input values $v_{i}$ and must be calculated at each predicate
call. This _dynamic_ part of the filter is the product of $k$ factors
$b_{1}b_{2}\ldots b_{k}$, where $b_{j}$ is the absolute value of either one of
the input variables $v_{i}$ or a difference $v_{a}-v_{b}$ between two of them.
Note that factors need not be different from each other (e.g. $b_{i}$ and
$b_{j}$ may be the same for $i\neq j$). Differently from [5], we slightly
overestimate $B(v_{1},...,v_{n})$ for the sake of efficient caching.
Specifically, we replace the product $b_{1}b_{2}\ldots b_{k}$ with the power
$\beta^{k}$, where $\beta=argmax\\{b_{1},\ldots,b_{k}\\}$. The reason for this
is the following.

Let $P(v_{1},\ldots,v_{k},i_{1},\ldots,i_{n})$ be a predicate where the
$v_{i}$’s are explicit parameters whereas the $i_{j}$’s are implicit. Any of
the implicit parameters $i$ is defined as a composition of explicit values
$w_{1}\ldots w_{m}$. Upon evaluation using the FP model, we first use the
$w$’s to calculate the $\lambda$’s and $d$’s, and then use these $\lambda$’s
and $d$’s together with the explicit parameters $v$’s to evaluate the whole
predicate. The semi-static filter for this composed calculation must take all
the variables into account. Thus, the dynamic component of our filter is
$\beta^{k}$, where $\beta$ is the maximum of the $b_{i}$’s deriving from
explicit variables $v$ and the $b_{j}$’s deriving from the $w$’s. The maximum
over the $b_{j}$’s is computed only once for each implicit point and cached
for later use. Hence, replacing the product with the power allows using a
single variable to cache the bounds for any $\lambda$ and $d$. When compared
with the caching benefits, the performance loss due to our overestimation is
virtually negligible. Fig. 3 shows an example of the code produced by our
tool.

Figure 3: Code produced by our tool for an indirect version of the orient2d
predicate. The first of the three points is implicitly defined by the
intersection of two lines that, in turn, are defined by explicit points
$w_{1},\ldots,w_{4}$. The values for $\lambda_{x}$, $\lambda_{y}$ and $d$
(l1x, l1y, d1) are calculated externally and passed to this function along
with the maximum $\beta$ deriving from the $w$’s (prev_maxval). The constant
value of $\delta(1)$ ( = 1.048458195263004e-13) accumulates the error bounds
of this function and those of the external function that calculates the
$\lambda$’s and $d$.


# Emotion-Driven Melody Harmonization
via Melodic Variation and Functional Representation

###### Abstract

Emotion-driven melody harmonization aims to generate diverse harmonies for a
single melody to convey desired emotions. Previous research found it hard to
alter the perceived emotional valence of lead sheets only by harmonizing the
same melody with different chords, which may be attributed to the constraints
imposed by the melody itself and the limitation of existing music
representation. In this paper, we propose a novel functional representation
for symbolic music. This new method takes musical keys into account,
recognizing their significant role in shaping music’s emotional character
through major-minor tonality. It also allows for melodic variation with
respect to keys and addresses the problem of data scarcity for better emotion
modeling. A Transformer is employed to harmonize key-adaptable melodies,
allowing for keys determined in rule-based or model-based manner. Experimental
results confirm the effectiveness of our new representation in generating key-
aware harmonies, with objective and subjective evaluations affirming the
potential of our approach to convey specific valence for versatile melody.

Index Terms—  Melody harmonization, symbolic music generation, data
representation, musical keys, emotion

##  1 Introduction

Inspired by the remarkable achievements made in symbolic music generation [1,
2, 3, 4], there is a growing interest in controlling high-level musical
features during the generation process, especially the perceived emotions from
music. Recent years have witnessed many efforts in unconditional music
generation [5, 6, 7, 8, 9] and melody harmonization [10, 11] to condition
their generation on emotion.

According to Russell’s famous Circumplex model of affect [12], emotion could
be represented in a two-dimensional space defined by valence and arousal,
where _valence_ related to the positiveness of an emotion and _arousal_ refers
to energy or activation [13]. Although many works are capable to control the
arousal of music, few of them succeeded in controlling the perceived valence.
For example, the piano music generation model in EMOPIA [5] fails in
generating low valence (i.e., negative) music, and the melody harmonization
model LHVAE [11] found it hard to change the overall emotion of music
(negative, neutral and positive) by only altering chords.

The ignorance of musical keys when modeling music data may be responsible for
the poor valence control. Valence is often found to be related to major-minor
tonality [14], and keys play important roles in affecting such tonality. The
histogram of keys derived from the emotion-labeled music dataset EMOPIA [5]
(Fig. 1) provides further support from a data perspective, where the
distribution skews to major keys for high valence clips and opposite trend for
low valence ones. However, to the best of our knowledge, none of existing
music generation works attempted to model keys explicitly.

Fig. 1: Key histogram of high/low valence clips from EMOPIA [5].

In this work, we focus on the _emotion-driven melody harmonization_ task,
which aims to convey desired emotions through harmonizing diverse chord
progressions for a single melody, resulting in the creation of lead sheets.
This simplified music format free from performance variances enable us to dive
into the impact of keys and chords on emotions, and we only consider the
valence aspect of emotion since arousal is usually related to performance-
level attributes and accompaniment patterns [15]. This task is more
challenging than conventional melody harmonization from three aspects.
Firstly, symbolic music datasets with both high-quality emotion and chord
labels are relatively scarce. Secondly, there are even less training examples
where one melody is accompanied with multiple chord progressions to convey
varying emotions. Thirdly, because the key is predetermined given a melody
line, there is limited room for chord progressions to contribute to the target
emotion. While the first challenge was partially addressed by previous works
[10, 11], the later two hinder the further improvement of valence control.

Observing the above, we propose a novel functional representation designed as
an alternative to REMI [1], a popular event representation that uses note
pitch values and chord names to encode symbolic music. Our method represents
both melody notes and chords with _Roman numerals_ relative to musical keys, a
_functional_ format considering the relationships between notes, chords and
scales (major or minor) [16]. In the harmonization process, driven by an
emotional condition (positive or negative), a key is determined in a rule-
based or model-based manner. Subsequently, a melody line encoded in functional
representation is accompanied to generate key-aware functional harmonies, a
process facilitated by a Transformer model.

Compared to note pitch values and chord names, which necessitate models to
infer keys implicitly, functional representation empowers us to explicitly
inform musical keys while modeling music pieces. Additionally, as melodies
across various scales are designed to be represented by the same set of
symbols for twelve scale degrees, the likelihood of encountering two melodies
with similar representations increases, and their accompanying chord
progressions and emotion labels can be regarded as paired samples, offering a
solution to the issue of data scarcity discussed above. Our representation
also supports transposition between parallel
keys111https://en.wikipedia.org/wiki/Parallel_key for any music pieces.
Specifically, when transposing to a parallel key, the scale degrees of notes
and chords remain consistent, while pitches may be adjusted to reflect the new
key mode, i.e., _melodic variation_. Through experiments, we aim to answer the
two research questions:

  * •

RQ#1: Can the proposed representation effectively model musical keys and yield
satisfactory harmonization outcomes?

  * •

RQ#2: Is it possible to generate different music variants from a single melody
to influence the perceived valence?

Our model implementations and checkpoints are open-source222Code:
https://github.com/Yuer867/EMO_Harmonizer, and generated samples could be
found in the demo webpage333Demo: https://yuer867.github.io/emo_harmonizer/.

Fig. 2: Illustration of (a) REMI [1] and (b) the proposed functional
representation, differing in note pitch and chord name events.

##  2 Related Work

Various model architectures and chord control methods have been explored in
melody harmonization [17, 18, 19, 20] and yielded remarkable results. In
contrast, progress in emotion-driven melody harmonization is limited due to
its inherent challenges. Some approaches relied on crowd-sourced emotion tags
as weak labels to obtain music samples with both emotion and chord labels
[10], while some approaches manually calculated valence simply from chord
labels [11]. For better emotion annotations, we use a pop piano performance
dataset with human-annotated emotions [5] and apply rule-based algorithm to
extract chord labels. We also leverage a high-quality lead sheet dataset for
pre-training. However, existing methods still struggle with controlling
emotions effectively, as the discrepancy between the perceived and target
emotions remain large. Moreover, none of these works consider the influence of
musical keys in harmonization. They typically transpose all samples to C major
/ c minor for simplicity [17, 18], or to all 12 keys for augmentation [20].
These challenges and constraints motivate us to study functional
representation for better emotion controllability.

##  3 Method

###  3.1 Functional Representation

The proposed functional representation is designed based on REMI [1], a widely
used event (token) based representation for symbolic music, but with different
note and chord events assisting to model the emotion and key information
better. See Fig. 2 for illustrations.

Emotion Events To denote distinct valences and affect overall properties, we
begin the event sequence with Emotion_Positive or Emotion_Negative to indicate
the emotional character of music clips, following the approach introduced in
CTRL [21].

Key Events After the Emotion_* event, a KEY_* event is appended to indicate
the key property. A total of 24 types (12 tonic notes with two modes each) are
used in this work.

Bar, Sub-Beat and EOS Events The same as REMI, a BAR event is used when a new
bar begins, a SUB-BEAT_* event points to one of 16 possible discrete locations
in a bar, and an EOS event will end the whole lead sheet.

Chord Events In existing melody harmonization studies, chord names are
commonly employed to denote root notes and chord qualities. For example, Fmaj
represents the chord F-A-C with root F and major quality. These symbols are
useful when all music pieces share the same tonic, i.e., C major and c minor,
but they overlook the variations in chord functions of the same chord across
different keys due to different scale degrees. For example, while Fmaj serves
tonic function in F major, it takes on subdominant function in C major.
Moreover, the chord progressions following specific functional harmony rules
establish tonality and convey musical emotion [16].

To introduce chord functions for better harmonization, we adopt Roman numerals
adopted from Roman Numeral Analysis [22] to notate chord roots independent of
keys (see Fig. 3). Given the key event, root notes belonging to the key scale
(such as C, D, E, F, G, A and B in C major) are directly converted into Roman
numerals based on their scale degrees relative to the tonic. For roots outside
the scale, we employ a direct conversion for I#, II#, IV#, V# and VI#
appearing in major keys, but randomly assign III# and VII# who only appear in
minor keys as one of their neighboring degrees. The conversion from Roman
numerals back to letters follows a similar process. This design ensures the
notation remaining key-independent and allows for harmonious transposition to
any desired key while minimizing pitch variations resulting from the
conversion. The notations of chord qualities remain unchanged. A CHORD_* event
appears every beat, even when there are no chord changes.

Fig. 3: Illustration of the conversion between letters and Roman numerals in
the cases of C major / c minor. The solid line represents a direct one-to-one
conversion, while the dotted line stands for a random conversion to either one
of them.

Note-related Events In REMI, every single note is denoted by Pitch_* and
Duration_* events. The former indicates the onset of pitches, ranging from
21(A0) to 108(C8). To enable key-aware pitches and melodic variation, we
decompose Pitch_* event into Octave_* and Degree_* events, where degrees are
computed as chord roots and octaves match the tonic’s octave in the key scale.
For example, pitch 63(D#4) will be decomposed into Octave_4 and Degree_III in
C minor. Note that this work omits performance-related events such as
VELOCITY_* and TEMPO_*.

Fig. 4 shows the functional representations of two bars taken from different
songs in EMOPIA. Two melody lines are represented in the same formats under
their corresponding keys (‘D’ major and ‘c’ minor), but harmonized with
different chords to convey distinct emotions (‘P’ositive and ‘N’egative). The
one on the right hand side can be considered as a melodic variation of the
left one.

Fig. 4: Examples of functional representations for two bars, featuring solely
emotion, key, note pitches, and chords for simplification.

###  3.2 Model Architecture

We follow the conditional generation framework proposed in Compound Word
Transformer [23] for our task. We firstly predict a key event $k$ conditioned
on an emotion event $e$, and then generate a chord sequence $C$ by
accompanying a melody sequence $M$ with $e$ and $k$, which could be formalized
as $p(k,C|e,M)=p(k|e)p(C|e,k,M)$. A sequence-to-sequence model [23, 24] is
applied to learn $p(k|e)$ and $p(C|e,k,M)$ simultaneously. Moreover, with the
positions of BAR events, $M$ and $C$ are further segmented into
$\\{M_{1},\cdots,M_{b}\\}$ and $\\{C_{1},\cdots,C_{b}\\}$, where $b$ is the
number of bars. The segmented sequences are interleaved in the form of
$\\{\cdots$TRACK_MELODY, $M_{i}$, TRACK_CHORD, $C_{i}\cdots\\}$ with
additional TRACK_* events, so that the target chord bar $C_{i}$ is nearest to
its corresponding conditions $M_{i}$ and the dependency between melody and
chord is easier to learn. The final model is summarized as
$p(k,C|e,M)=p(k|e)\prod_{i=1}^{b}p(C_{i}|e,k,M_{\leq i},C_{<i})$ and minimized
the negative log-likelihood loss of the generated sequence.

At inference time, given an emotion condition $e$, a key event $k$ could be
determined in rule-based or model-based manner. In the former approach, we
enforce the original key to its parallel major key for positive emotion or
parallel minor key otherwise, inspired by Fig. 1. However, a song in major key
could also convey negative emotion, so alternatively, the key can be predicted
from $p(k|e)$ learned above. Note that the corresponding melody will be
adjusted as key changes.

As the size of the dataset with emotion labels is not big enough, we pretrain
the model with a large lead sheet dataset without emotion annotations to
establish a robust understanding of relationship between melody and chord,
where EMOTION_NONE is used as the emotion event. We then finetune the model on
EMOPIA [5] to learn harmonization styles specific to different emotion
conditions.

##  4 Experiments

###  4.1 Datasets, Preprocessing, and Model Settings

We adopt a large lead sheet dataset collected from HookTheory [25] released in
SheetSage [26] for pre-training (refered to as “HookTheory” hereafter). Each
piece in HookTheory contains high-quality, human-transcribed melody, chord and
key annotations. After removing the pieces that are not in 4/4 time signature
or major / minor keys, over 18k segments remain. We simplify 249 chord quality
classes to 11 types including major, minor, augment, diminish, suspend2,
suspend4, major7, minor7, dominant7, diminish7, half-diminish7, which may be
sufficient to convey most emotional characters.

We use the piano MIDI dataset EMOPIA for model finetuning. It contains 1,071
music clips with human-annotated emotion labels. A four-class taxonomy adopted
from Russell’s model [12] is used for annotations, including HVHA (high
valence high arousal), HVLA (high valence low arousal), LVHA (low valence high
arousal) and LVLA (low valence low arousal). Since we only consider the
valence dimension, clips in HVHA and HVLA are combined into Positive class and
others into Negative class. To form the lead sheet samples from piano
performance, we adopt a method similar to the one proposed in [27] to extract
chord labels in 11 quality types, and apply a heuristic rule-based method to
extract the melody line [28].

The statistics of two datasets after preprocessing are shown in Table 1. The
clips in HookThoery are randomly divided into train and validation splits, and
we follow the stratified split code provided by EMOPIA for train and
validation, both with the ratio of 9:1. In our data representations, the
vocabulary size of events is 217.

Dataset | # clips | # bars (avg.) | # events (avg.)
---|---|---|---
HookTheory [26] | 18,206 | 10.84 | 339.27
EMOPIA [5] | 1,071 | 16.96 | 591.80

Table 1: Summary of datasets used in our experiments.

A 12-layer linear Transformer with Performer attention [29] is used for
generation (8 attention head, 512 hidden state dim., 38 million parameters)
and trained with a batch size of 4 and a maximum sequence length of 1,024
(longer than 99.2% and 94.6% of clips in HookTheory and EMOPIA respectively)
on a Tesla V100 GPU with 32G VRAM. We use Adam optimizer with 200 steps of
warmup to maximum learning rate of 1e–4 and 1e–5 for pretrain and finetune
respectively, both followed by 500k steps of cosine decay. Checkpoints with
lowest validation loss are used for inference, using nucleus sampling [30]
with temperatured softmax ($\tau=1.1$, $p=0.99$).

###  4.2 Objective Evaluation

We firstly study RQ#1, i.e., can the proposed representation effectively model
musical keys and yield satisfactory harmonization outcomes. We consider two
baselines for comparison. One represents music clips in REMI [1] directly
without any other processing, following the approach used in emotion-
conditioned melody harmonization work LHVAE [11]. The other transposes
(‘trans’) all clips to C major / c minor before REMI encoding, standing for a
common method of existing melody harmonization works [17, 18]. We do not
compare with other model architectures as our primary focus is to study the
impact of different representations. For functional representation, besides
the full format, we also compare one ablated version (‘ablated’), where chords
are encoded in the functional format but melody notes are represented by note
pitch values as in REMI.

For each of the variants, we harmonize all melodies in validation set under
the conditions of positive and negative emotion respectively, i.e., $88\times
2$ samples each variant, to check their harmonization performance. Two groups
of objective evaluation metrics are considered, including three metrics
proposed in [17] to evaluate the harmonicity between chord and melody: chord
tone to non-chord tone ratio (CTnCTR), pitch consonance score (PCS), melody-
chord tonal distance (MCTD), and two metrics _newly_ proposed here for
evaluating the ability of modeling keys:

  * •

Root ratio (RR): the ratio of chord _roots_ in key scale.

  * •

Note ratio (NR): the ratio of chord _notes_ in key scale.

Table 2 shows that, while encoding music clips into REMI representation
directly have poor performances in all metrics, providing musical keys
information by transposing to C major / c minor or adding key events with
functional representation could greatly improve their harmonicity and almost
meet the real data. This indicates that musical key is critical to melody
harmonization, and that our representations enables efficient learning of
their relationships despite of fewer training data for each key type.
Moreover, the ablated version suffer losses on melody-related metrics, showing
the necessity of encoding melodies in functional formats. In short, the answer
to RQ#1 is that when applying functional representation on both melody notes
and chord labels, musical keys could be effectively modeled and the melodies
are well harmonized.

Methods | CTnCTR | PCS | MCTD | RR | NR
---|---|---|---|---|---
REMI [1, 11] | 0.291 | 0.211 | 1.731 | .647 | .592
REMI (trans) | 0.747 | 1.557 | 1.347 | .927 | .923
Ours | 0.750 | 1.487 | 1.343 | .943 | .935
Ours (ablated) | 0.530 | 0.880 | 1.501 | .952 | .942
Real data | 0.801 | 1.613 | 1.314 | .935 | .926

Table 2: Melody harmonization objective evaluation results (the closer to the
real data, the better).

Methods | Objective $\downarrow$ | Subjective $\uparrow$
---|---|---
QD | PD | Novice | Expert | Total
REMI (trans) | .0128 | .0218 |  $-$0.04$\pm$0.92 |  $-$0.33$\pm$0.94 |  $-$0.16$\pm$0.94
REMI (rule) | .0140 | .0161 |  0.44$\pm$1.12 |  0.69$\pm$1.08 |  0.54$\pm$1.11
Ours | .0127 | .0119 | — | — | —
Ours (rule) | .0084 | .0116 |  0.61$\pm$0.91 |  0.58$\pm$1.11 | 0.60$\pm$1.00
Ours (model) | .0114 | .0130 |  $-$0.093$\pm$1.36 |  0.39$\pm$1.16 |  0.10$\pm$1.31

Table 3: Emotion controllability results via objective evaluation (the lower
the better) and subjective evaluation (the higher the better).

Next, we study RQ#2, i.e., whether we could control the valence when
generating different music variants from a melody. We select the methods
perform well above, i.e., REMI (trans) and our functional representation.
Besides simply setting the emotion event to the target emotion without
changing keys, we further examine rule-based (‘rule’) and model-based
(‘model’) methods discussed in Section 3.2 to determine keys conditioned on
emotions, and build five models in total. Note that the ‘rule’ variant of REMI
(trans) is implemented by transposing the melody to major / minor keys
directly. To quantify how well the generated samples conform to the emotion
conditions, we propose the following two metrics:

  * •

Quality Distribution (QD): compute the KL divergence between the chord quality
distributions of generated samples and real data for positive and negative
emotions respectively, then take average.

  * •

Progression Distribution (PD): calculate the KL divergence between the bi-gram
chord progression distributions as above and compute the average. Progression
refers to the difference between two consecutive chord root in chromatic
scale. For example, the progression between Dmaj and Fmin is 3.

To obtain reliable distributions, for each combination, we generate harmonies
for all melodies in validation set five times under both positive and negative
emotion conditions, i.e., $88\times 5\times 2$ samples, to match the size of
training samples. As is shown in the left part of Table 3, the combination of
functional representation and rule-based key determination achieves the lowest
distribution distances for two metrics. Furthermore, all the versions using
functional representation performs better than REMI ones, indicating the
effectiveness of this representation to model emotion-related properties of
music.

###  4.3 User Study

An online survey was deployed to collect user responses. Every subject needs
to listen to 16 music pieces, including two original pieces in positive and
two in negative as well as their four variants generated with the same melody
but the opposite emotion condition by four methods444The method using
functional representation without key changes is removed here to lower the
burden of users introduced above. These original pieces are randomly drawn
from the validation set. For each group of samples, users will assess how the
variants (in random order) compare to the original piece in terms of their
emotional differences on a five-point scale, specifically whether they convey
a much more positive(2), more positive(1), unchanged(0), more negative(–1) or
much more negative(–2) emotion, without knowing the pre-defined emotion
conditions. This design is inspired by the finding that it is easier for human
subjects to make relative valence comparison rather than to assign absolute
rating [31]. 23 subjects participated in the survey (‘Total’), 9 of them with
$\geq$4 years of musical training or experience (‘Expert’).

The right side of Table 3 shows the average scores from the subjects, with a
higher score ($\in[-2,2]$) indicating better emotion controllability, i.e.,
the perceived valence matches the given condition. Totally speaking, the
combination of functional representation and rule-based key determination
performs best, followed by REMI with rule-based keys, which indicates that
musical keys play a significant role in influencing the perceived valence
through melodic variation and the functional representation generates slightly
better harmonies to support the perception of desire emotions, which answers
the RQ#2. Moreover, users are still unable to perceive the emotion conversions
by simply harmonizing the original melody without any key changes, aligning
with the previous findings [11]. When determining keys in model-based method,
sometimes major (minor) keys will be sampled for negative (positive) emotions
and the generated harmonies lack the ability to obviously change the perceived
valence by themselves, yielding worse emotion controllability. Namely, the
rule-based way seems to outperform the model-based way in determining keys.
Improving the latter will be a focus of future work.

Fig. 5: Comparison of subjective results under different settings.

When diving into user responses from different musical backgrounds, those with
longer years of music-related experiences seem to favor the REMI-generated
samples, while people with less experiences choose the functional
representation. Intuitively, with more musical training, people may analyze
the conveyed emotions more from the perspective of music theory, while novices
simply rely on their emotional senses, but it is hard to say which method is
more accurate. Moreover, if examining the cases of ‘negative to positive’
(i.e., generating positive variants for negative original pieces) and
‘positive to negative’ separately (Fig. 5), the control of low valence is
better with more samples match the negative condition, while half of samples
conditioned on positive are still identified as negative ones as their ground-
truth emotions. It seems that minor keys have significant influence on
negative feelings, while the melody flows have strong influences on positive
moods.

##  5 Conclusion and Future Work

In this paper, we have proposed a novel functional representation for symbolic
music, which represents melody notes and chords with Roman numerals relative
to musical keys. A Transformer-based framework is then adopted to harmonize
melodies conditioned on emotional valence. Objective assessments validate our
approach’s effectiveness in key modeling, while subjective evaluations confirm
its ability to convey desired emotional valence. Future endeavors may focus on
the improvement of model-based approach and the control of arousal aspect of
emotion through accompaniment generation.

## References

  * [1] Yu-Siang Huang and Yi-Hsuan Yang,  “Pop Music Transformer: Beat-based modeling and generation of expressive pop piano compositions,”  in Proc. ACM Multimedia, 2020.
  * [2] Cheng-Zhi Anna Huang, Ashish Vaswani, Jakob Uszkoreit, Ian Simon, Curtis Hawthorne, Noam Shazeer, Andrew M. Dai, Matthew D. Hoffman, Monica Dinculescu, and Douglas Eck,  “Music transformer: Generating music with long-term structure,”  in Proc. ICLR, 2019.
  * [3] Dimitri von Rütte, Luca Biggio, Yannic Kilcher, and Thomas Hofmann,  “FIGARO: Generating symbolic music with fine-grained artistic control,”  in Proc. ICLR, 2023.
  * [4] Peiling Lu, Xin Xu, Chenfei Kang, Botao Yu, Chengyi Xing, Xu Tan, and Jiang Bian,  “MuseCoco: Generating symbolic music from text,”  CoRR, vol. abs/2306.00110, 2023.
  * [5] Hsiao-Tzu Hung, Joann Ching, Seungheon Doh, Nabin Kim, Juhan Nam, and Yi-Hsuan Yang,  “EMOPIA: A multi-modal pop piano dataset for emotion recognition and emotion-based music generation,”  in Proc. ISMIR, 2021.
  * [6] Chenfei Kang, Peiling Lu, Botao Yu, Xu Tan, Wei Ye, Shikun Zhang, and Jiang Bian,  “EmoGen: Eliminating subjective bias in emotional music generation,”  CoRR, vol. abs/2307.01229, 2023.
  * [7] Shulei Ji and Xinyu Yang,  “EmoMusicTV: Emotion-conditioned symbolic music generation with hierarchical Transformer VAE,”  IEEE Transactions on Multimedia, 2023.
  * [8] Lucas Ferreira and Jim Whitehead,  “Learning to generate music with sentiment,”  in Proc. ISMIR, 2019.
  * [9] Dimos Makris, Kat R. Agres, and Dorien Herremans,  “Generating lead sheets with affect: A novel conditional seq2seq framework,”  in Proc. IJCNN, 2021.
  * [10] Takuya Takahashi and Mathieu Barthet,  “Emotion-driven harmonisation and tempo arrangement of melodies using transfer learning,”  in Proc. ISMIR, 2022.
  * [11] Shulei Ji and Xinyu Yang,  “Emotion-conditioned melody harmonization with hierarchical variational autoencoder,”  CoRR, vol. abs/2306.03718, 2023.
  * [12] James A. Russell,  “A circumplex model of affect,”  Journal of Personality and Social Psychology, 1980.
  * [13] Juan Sebastián Gómez Cañón, Estefanía Cano, Tuomas Eerola, Perfecto Herrera, Xiao Hu, Yi-Hsuan Yang, and Emilia Gómez,  “Music emotion recognition: Toward new, robust standards in personalized and context-sensitive applications,”  IEEE Signal Processing Magazine, 2021.
  * [14] Renato Panda, Ricardo Malheiro, and Rui Pedro Paiva,  “Audio features for music emotion recognition: A survey,”  IEEE Trans. Affective Computing, 2020.
  * [15] Yi-Chan Wu and Homer H. Chen,  “Generation of affective accompaniment in accordance with emotion flow,”  IEEE Trans. Audio Speech and Language Processing, 2016.
  * [16] Tsung-Ping Chen and Li Su,  “Functional harmony recognition of symbolic music data with multi-task recurrent neural networks,”  in Proc. ISMIR, 2018.
  * [17] Yin-Cheng Yeh, Wen-Yi Hsiao, Satoru Fukayama, Tetsuro Kitahara, Benjamin Genchel, Hao-Min Liu, Hao-Wen Dong, Yian Chen, Terence Leong, and Yi-Hsuan Yang,  “Automatic melody harmonization with triad chords: A comparative study,”  CoRR, vol. abs/2001.02360, 2020.
  * [18] Chung-En Sun, Yi-Wei Chen, Hung-Shin Lee, Yen-Hsing Chen, and Hsin-Min Wang,  “Melody harmonization using orderless nade, chord balancing, and blocked gibbs sampling,”  Proc. ICASSP, 2020.
  * [19] Yi-Wei Chen, Hung-Shin Lee, Yen-Hsing Chen, and Hsin-Min Wang,  “SurpriseNet: Melody harmonization conditioning on user-controlled surprise contours,”  in Proc. ISMIR, 2021.
  * [20] Jingwei Zhao, Gus Xia, and Ye Wang,  “Domain adversarial training on conditional variational auto-encoder for controllable music generation,”  in Proc. ISMIR, 2022.
  * [21] Nitish Shirish Keskar, Bryan McCann, Lav R. Varshney, Caiming Xiong, and Richard Socher,  “CTRL: A conditional transformer language model for controllable generation,”  CoRR, vol. abs/1909.05858, 2019.
  * [22] Gianluca Micchi, Mark Gotham, and Mathieu Giraud,  “Not all roads lead to Rome: Pitch representation and model architecture for automatic harmonic analysis,”  TISMIR, 2020.
  * [23] Wen-Yi Hsiao, Jen-Yu Liu, Yin-Cheng Yeh, and Yi-Hsuan Yang,  “Compound Word Transformer: Learning to compose full-song music over dynamic directed hypergraphs,”  in Proc. AAAI, 2021.
  * [24] Shih-Lun Wu and Yi-Hsuan Yang,  “Compose & Embellish: Well-structured piano performance generation via a two-stage approach,”  in Proc. ICASSP, 2023.
  * [25] “HookTheory,” https://www.hooktheory.com/ [Accessed: (September 1, 2023)].
  * [26] Chris Donahue, John Thickstun, and Percy Liang,  “Melody transcription via generative pre-training,”  in Proc. ISMIR, 2022.
  * [27] Bryan Pardo and William P. Birmingham,  “Algorithms for chordal analysis,”  Computer Music Journal, 2002.
  * [28] “Midi_Toolkit,” https://github.com/RetroCirce/Midi_Toolkit [Accessed: (September 1, 2023)].
  * [29] Krzysztof Marcin Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tamás Sarlós, Peter Hawkins, Jared Quincy Davis, Afroz Mohiuddin, Lukasz Kaiser, David Benjamin Belanger, Lucy J. Colwell, and Adrian Weller,  “Rethinking attention with Performers,”  in Proc. ICLR, 2021.
  * [30] Ari Holtzman, Jan Buys, Maxwell Forbes, and Yejin Choi,  “The curious case of neural text degeneration,”  in Proc. ICLR, 2019.
  * [31] Yi-Hsuan Yang and Homer H. Chen,  “Ranking-based emotion recognition for music organization and retrieval,”  IEEE Trans. Speech Audio Process, 2011.