# Efficient Computation of Topological Integral Transforms

Vadim Lebovici, Steve Oudot, Hugo Passe

► To cite this version:

Vadim Lebovici, Steve Oudot, Hugo Passe. Efficient Computation of Topological Integral Transforms. SEA 2024 - 22nd International Symposium on Experimental Algorithms, Jul 2024, Vienna, Austria. 10.4230/LIPIcs.SEA.2024.22 . hal-04856435

HAL Id: hal-04856435

https://inria.hal.science/hal-04856435v1

Submitted on 26 Feb 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient Computation of Topological Integral Transforms

Vadim Lebovici<fnref n="1" /> ☐

Mathematical Institute, University of Oxford, United Kingdom

Steve Oudot ☐

Inria Saclay, Palaiseau, France

Hugo Passe ☐

Ecole Normale Supérieure de Lyon, France

## Abstract

Topological integral transforms have found many applications in shape analysis, from prediction of clinical outcomes in brain cancer to analysis of barley seeds. Using Euler characteristic as a measure, these objects record rich geometric information on weighted polytopal complexes. While some implementations exist, they only enable discretized representations of the transforms, and they do not handle weighted complexes (such as for instance images). Moreover, recent *hybrid transforms* lack an implementation.

In this paper, we introduce eucalc, a novel implementation of three topological integral transforms – the Euler characteristic transform, the Radon transform, and hybrid transforms – for weighted cubical complexes. Leveraging piecewise linear Morse theory and Euler calculus, the algorithms significantly reduce computational complexity by focusing on critical points. Our software provides exact representations of transforms, handles both binary and grayscale images, and supports multi-core processing. It is publicly available as a C++ library with a Python wrapper. We present mathematical foundations, implementation details, and experimental evaluations, demonstrating eucalc's efficiency.

**2012 ACM Subject Classification** Mathematics of computing → Algebraic topology

**Keywords and phrases** Topological data analysis, Euler calculus, Topological integral transform, Euler characteristic transform, Hybrid transforms

**Digital Object Identifier** 10.4230/LIPIcs.SEA.2024.22

**Supplementary Material** Software (Source Code): https://github.com/HugoPasse/Eucalc
archived at swh:1:dir:cee8afe463352094a26fc586f6cc1c92913de409

**Funding** Vadim Lebovici : Partly funded by EPSRC EP/R018472/1.

# 1 Introduction

## 1.1 Motivations

Since its introduction in the late 80s by Viro [17] and Schapira [13], *Euler calculus* – the integral calculus with respect to the Euler characteristic – has been of increasing interest in topological data analysis. Using this integral calculus, one can define integral transforms conveying significant topological information. The most famous one, the *Euler characteristic transform* (ECT), has found many applications in shape analysis, for instance in the prediction of clinical outcomes in brain cancer [5], in the analysis of barley seeds [2] or in the recovery of morphological variations across genera of primates [15]. More formally, given a sufficiently

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹ For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.</footnote>
## 22:2 Efficient Computation of Topological Integral Transforms

tame compact shape $K$ in Euclidean space (e.g., an embedded polytopal complex), the ECT records the Euler characteristic of the intersection $\xi^{-1}(-\infty, t] \cap K$ for all linear forms $\xi: \mathbb{R}^n \to \mathbb{R}$ and all $t \in \mathbb{R}$. What is remarkable is that this transform entirely characterizes the shape $K$, that is, it is injective: if two shapes have same ECT, then they are equal [7, 8, 14, 16]. The ECT is an instance of a more general transform called the *Radon transform* [14] and the injectivity of the ECT is a consequence of Theorem 3.1 in loc. cit.. Another natural specialization of the Radon transform, simply called *Radon transform* in this paper, records the Euler characteristic of $\xi^{-1}(t) \cap K$ for all linear forms $\xi: \mathbb{R}^n \to \mathbb{R}$ and all $t \in \mathbb{R}$.

The most problematic aspect of Euler calculus is its instability under numerical approximations [6]. To remedy it, integral transforms mixing Euler calculus and Lebesgue integration were introduced in [9], then generalized under the name *hybrid transforms* (HT) in [11]. These transforms integrate the Radon transform against a chosen kernel to smooth out and compress the information. The use of a kernel provides a wide variety of potential summaries that can be used to emphasize different aspects of data [10]. All these transforms are defined more generally for any constructible function, but in this article we restrict ourselves to the case of weighted polytopal complexes, i.e., polytopal complexes for which each cell is assigned an integer.

## 1.2 Existing implementations and limitations

In this paper, we are interested in the computation of topological integral transforms. Several implementations already exist. The Python package `demeter` [1, 2] implements a discretized version of the ECT for axis-aligned cubical complexes. This package turns an $n$-dimensional binary image into an embedded cubical complex $C$ in $\mathbb{R}^n$ in a prescribed way: vertices are defined by non-zero pixels and higher dimensional cubes are defined by adjacency relations. The ECT of $C$ is then output for any direction $\xi: \mathbb{R}^n \to \mathbb{R}$ as a discretization of the *Euler characteristic curve* $t \in \mathbb{R} \mapsto \chi(\xi^{-1}(-\infty, t] \cap C)$ for some chosen discretization step $T$. The computation of the curve is done using the classical formula for the Euler characteristic $\chi(\xi^{-1}(-\infty, t] \cap C)$ as a sum of terms $(-1)^{\dim(C)}$ over all cells $c$ of $C$ such that $c \subseteq \xi^{-1}(-\infty, t]$. A bucket sort of $\xi$ values on the vertices optimizes the computation of the discretization. The Python package `Sinatra-Pro` [15] implements a similar algorithm for triangular meshes both for the ECT and its smoothed version (SECT) [5].

These implementations suffer from several essential limitations. First of all, they only work with polytopal complexes, not weighted polytopal complexes. This prevents, e.g., the processing of images, which are ubiquitous in applications. Furthermore, these implementations only output discretized Euler characteristic curves. This is important because, while $2^n$ well-chosen Euler characteristic curves are sufficient to fully reconstruct a cubical complex embedded in $\mathbb{R}^n$, each curve must be computed exactly for the reconstruction to be correct. Reconstructing the complex from discretized curves will likely lead to substantial errors in the reconstruction, due to the instability of Euler calculus with respect to numerical approximations. On the contrary, an exact computation of integral transforms would pave the way for the implementation of their left-inverses. A software computing an exact version of the ECT on weighted cubical complexes together with an algorithm for its left-inverse does exist, but only for two- and three-dimensional complexes [4]. Besides, this software (as all aforementioned ones) does not allow for the computation of the Radon and hybrid transforms. There is only one recently released software [10] that computes hybrid transforms for some specific complexes arising in the context of multi-parameter persistence. However, no algorithm exists for general polytopal complexes, nor even for cubical ones. To the best of our knowledge, no implementation of the Radon transform exists.
## 1.3 Contributions

We introduce fast implementations of the ECT, Radon transform, and hybrid transforms for weighted (axis-aligned) cubical complexes built from grayscale images. After a pre-processing step, our implementation can evaluate the transform exactly at any vector $\xi$ in time at most linear (and often largely sublinear in practice) in the number of vertices in the complex. This is thanks to two major optimizations.

Our first and main optimization comes from the use of piecewise linear Morse theory [3], which guarantees that, under some assumptions that are automatically satisfied by linear forms $\xi: \mathbb{R}^n \to \mathbb{R}$ on polytopal complexes $K$, the topology of $\xi^{-1}(-\infty, t] \cap K$ changes only when $t \in \mathbb{R}$ belongs to a finite set of so-called *critical values*. These critical values only appear when $\xi^{-1}(t)$ contains a vertex of the complex called a *critical point*. The set of critical values of $\xi$ together with the associated changes in the Euler characteristic of sublevel-sets yield a finite and exact representation of the associated Euler characteristic curve – see for instance [7] for a precise statement in the case of unweighted cubical complexes. This reduces the computation from a sum over all cells to a sum over only a few critical points. As an illustration, we show that there are between 34 and 21 times fewer critical points than vertices on the images of the *fashion_MNIST* data set.

Our second optimization comes from the observation that two linear forms $\xi$ and $\xi'$ whose canonical coordinates have same signs share the same critical points on axis-aligned cubical complexes. Precomputing these points for $2^n$ well-chosen directions enables fast evaluation of topological transforms in many directions at a very low cost. Moreover, the formalism of Euler calculus suggests that the results of piecewise linear Morse theory also hold for weighted cubical complexes. We exploit this remark in our implementation. Consequently, computing topological transforms of grayscale images is done in less than twice the time for binary ones.

Our software, called *eucalc* [12], is available as a C++ library together with a Python wrapper to make it easily compatible with scikit-learn like interfaces. It provides a data structure for weighted embedded cubical complexes allowing for fast preprocessing of critical points. Our data structure builds on the implementation of cubical complexes from the *GUDHI* library to benefit from their flexibility of instanciation. Moreover, our structure represents exact versions of the ECT and the Radon transform, together with evaluation and discretization routines. Our implementation also includes a parallelized version of the preprocessing step and the computation of the hybrid transform.

# 2 Preliminaries

Let $n \ge 2$ be an integer. The space of linear forms on $\mathbb{R}^n$ will be canonically identified with $\mathbb{R}^n$. For clarity, we refer to linear forms with greek letters $\xi \in \mathbb{R}^n$, and we denote the canonical dot product of $\xi$ and $x \in \mathbb{R}^n$ by $\langle \xi; x \rangle$. For any function $f: \mathbb{R}^n \to \mathbb{R}$ and $t \in \mathbb{R}$, we denote by $\{f = t\}$ and $\{f \le t\}$ the level set $\{x \in \mathbb{R}^n : f(x) = t\}$ and the sublevel set $\{x \in \mathbb{R}^n : f(x) \le t\}$ respectively.

## 2.1 Weighted polytopal complexes

While our algorithms deal with cubical complexes, we introduce our main objects in the more general setting of polytopal complexes for the sake of clarity. We refer to Ziegler [18] for more details on polytopal complexes. We call *polytope* a bounded intersection of a finite number of closed half-spaces of $\mathbb{R}^n$. Equivalently, a polytope is the convex hull of a finite
22:4 Efficient Computation of Topological Integral Transforms

number of points in $\mathbb{R}^n$, called its *vertices* [18, Thm. 1.1]. We denote by Vert($P$) the (finite) set of all vertices of a polytope $P$, and we call it its *vertex set*. A *face* of $P$ is any set of the form $F = P \cap \{\xi = t\}$ where $\xi \in \mathbb{R}^n$ and $t \in \mathbb{R}$ are such that $P \subseteq \{\xi \le t\}$. The *dimension* of a polytope is the dimension of the smallest affine subspace of $\mathbb{R}^n$ containing it. A *polytopal complex* $\mathcal{P}$ is a finite collection of polytopes in $\mathbb{R}^n$ (called *cells*) such that (i) the empty cell is in $\mathcal{P}$, (ii) if $P$ is in $\mathcal{P}$, then all the faces of $P$ are also in $\mathcal{P}$, and (iii) if $P$ and $Q$ are in $\mathcal{P}$, then their intersection is a face of both $P$ and $Q$. One can always extend a collection of polytopes satisfying (i) and (iii) to a polytopal complex by adding all faces of the polytopes in the collection. This polytopal complex is the smallest one containing the polytopes of the given collection.

In this paper, we are mainly interested in cubical complexes. A *cube* in $\mathbb{R}^n$ is a subset $P = [a_1, b_1] \times \dots \times [a_n, b_n] \subseteq \mathbb{R}^n$ such that $b_i - a_i = b_j - a_j$ if $a_i \neq b_i$ and $a_j \neq b_j$. A *cubical complex* $\mathcal{P}$ is a polytopal complex whose polytopes are cubes in $\mathbb{R}^n$.

Let $\mathcal{P}$ be a polytopal complex in $\mathbb{R}^n$ and $\xi \in \mathbb{R}^n$. The *star* of a polytope $P \in \mathcal{P}$ is the collection $\text{St}(P)$ of polytopes $Q \in \mathcal{P}$ such that $Q \supseteq P$. For any $x \in \text{Vert}(\mathcal{P})$, the *lower star* of $x$ with respect to $\xi$ is the collection $\text{LwSt}(\xi, x)$ of polytopes $Q \in \text{St}(x)$ such that $\max_Q \xi = \xi(x)$. Similarly, we define the *upper star* $\text{UpSt}(\xi, x)$. These collections of polytopes are not complexes as they are not closed under taking faces.

Let $\mathcal{X}$ be a finite collection of polytopes in $\mathbb{R}^n$ and let $\varphi : \mathcal{X} \to \mathbb{Z}$. If $\mathcal{X}$ is a polytopal complex, then we call $\varphi$ a *weighted polytopal complex*. The *Euler characteristic* of $\varphi$ is defined as $\chi(\varphi) = \sum_{P \in \mathcal{X}} (-1)^{\dim(P)} \varphi(P)$. This definition generalizes the usual notion of Euler characteristic of a polytopal complex [18, Corollary 8.17] and is a particular case of Euler integration [17, 13].

## 2.2 Integral transforms

Let $\xi \in \mathbb{R}^n$ and $t \in \mathbb{R}$. Consider the collections of polytopes $\{P \cap \{\xi \le t\} : P \in \mathcal{P}\}$ and $\{P \cap \{\xi = t\} : P \in \mathcal{P}\}$ and denote respectively by $\mathcal{P}_{\{\xi \le t\}}$ and $\mathcal{P}_{\{\xi = t\}}$ the smallest polytopal complexes containing these collections. Any function $\varphi : \mathcal{P} \to \mathbb{Z}$ induces a function $\varphi_{\{\xi \le t\}} : \mathcal{P}_{\{\xi \le t\}} \to \mathbb{Z}$ defined for any $Q \in \mathcal{P}_{\{\xi \le t\}}$ by $\varphi_{\{\xi \le t\}}(Q) = \varphi(P)$ where $P$ is the smallest polytope of $\mathcal{P}$ containing $Q$. In that case, one has $\dim(Q) = \dim(P) - 1$. Similarly, the function $\varphi$ induces a function $\varphi_{\{\xi = t\}} : \mathcal{P}_{\{\xi = t\}} \to \mathbb{Z}$.

Let $\varphi : \mathcal{P} \to \mathbb{Z}$. The first topological integral transform we consider in this paper is the *Radon transform* [14], namely, the function $\mathcal{R}_{\varphi} : \mathbb{R}^n \times \mathbb{R} \to \mathbb{Z}$ defined for any $(\xi, t) \in \mathbb{R}^n \times \mathbb{R}$ by:

$$ \mathcal{R}_{\varphi}(\xi, t) = \chi(\varphi_{\{\xi=t\}}). $$

The second integral transform, the *Euler characteristic transform* [16], is the function $\text{ECT}_{\varphi} : \mathbb{R}^n \times \mathbb{R} \to \mathbb{Z}$ defined for any $(\xi, t) \in \mathbb{R}^n \times \mathbb{R}$ by:

$$ \text{ECT}_{\varphi}(\xi, t) = \chi(\varphi_{\{\xi \le t\}}). $$

Following [11], we consider specific smoothings of topological transforms. Let $\kappa : \mathbb{R} \to \mathbb{C}$ be locally integrable. The *hybrid transform* of $\varphi$ with kernel $\kappa$ is the function $\text{HT}_{\varphi} : \mathbb{R}^n \to \mathbb{C}$ defined for any $\xi \in \mathbb{R}^n$ by:

$$ \text{HT}_{\varphi}(\xi) = \int_{\mathbb{R}} \kappa(t) \mathcal{R}_{\varphi}(\xi, t) dt. $$
## 2.3 Naive algorithm

The following lemma is an easy consequence of the definitions.

**Lemma 1.** Let $\bar{\kappa} : \mathbb{R} \to \mathbb{C}$ be a primitive of $\kappa$.

$$ \mathcal{R}_{\varphi}(\xi, \cdot) = \sum_{P \in \mathcal{P}} \varphi(P) \mathbf{1}_{[\min_P \xi, \max_P \xi]}, $$

$$ \mathrm{ECT}_{\varphi}(\xi, \cdot) = \sum_{P \in \mathcal{P}} \varphi(P) \mathbf{1}_{[\min_P \xi, +\infty)}. $$

As a consequence,

$$ \mathrm{HT}_{\varphi}(\xi) = \sum_{P \in \mathcal{P}} \varphi(P) \cdot (\bar{\kappa}\left(\max_P \xi\right) - \bar{\kappa}\left(\min_P \xi\right)). $$

This lemma straightforwardly suggests naive algorithms to compute topological integral transforms over any polytopal complex. Doing so, the time complexity of computing any one of the above integral transforms over $\mathcal{P}$ is $O(N \cdot \#\mathcal{P})$ where $N$ is the maximum number of vertices of cells in $\mathcal{P}$. These naive algorithms will serve as a baseline for our future optimizations. Notice that $\#\mathcal{P}$ is exponential in the dimension.

## 3 Euler critical points

In this section, we slightly generalize the notion of Euler critical value of polytopal complexes [7, Def. 6.2] to weighted polytopal complexes. All mathematical results are folklore. We state them and provide proofs in Appendix B for completeness.

Let $\mathcal{P}$ be a polytopal complex, $\xi \in \mathbb{R}^n$ and $\varphi : \mathcal{P} \to \mathbb{Z}$. Further assume that $\xi$ is generic on $\mathcal{P}$, i.e., no polytope of $\mathcal{P}$ of positive dimension is contained in a level set of $\xi$. We define the ordinary and classical critical values of $\xi$ at a vertex $x$ of $\mathcal{P}$ respectively by:

$$ \Delta_{\varphi}^{\mathrm{ord}}(\xi, x) = \chi(\varphi|_{\mathrm{UpSt}(\xi, x)}) - \chi(\varphi|_{\mathrm{LwSt}(\xi, x)}), $$

$$ \Delta_{\varphi}^{-}(\xi, x) = \chi(\varphi|_{\mathrm{LwSt}(\xi, x)}). $$

A vertex $x$ of $\mathcal{P}$ is an Euler $\varphi$-critical point of $\xi$ if one of $\Delta_{\varphi}^{-}(\xi, x)$ or $\Delta_{\varphi}^{\mathrm{ord}}(\xi, x)$ is non-zero. Such a point is called ordinary when $\Delta_{\varphi}^{\mathrm{ord}}(\xi, x) \neq 0$ and classical when $\Delta_{\varphi}^{-}(\xi, x) \neq 0$. Euler $\varphi$-critical points are sometimes called critical points when $\varphi$ and $\xi$ are clear from the context. Note that a critical point can be both classical and ordinary. The set of classical critical points is denoted by $\mathrm{Crit}_{\varphi}^{-}(\xi)$, the set of ordinary critical points by $\mathrm{Crit}_{\varphi}^{\mathrm{ord}}(\xi)$ and the set of all critical points by $\mathrm{Crit}_{\varphi}(\xi) = \mathrm{Crit}_{\varphi}^{-}(\xi) \cup \mathrm{Crit}_{\varphi}^{\mathrm{ord}}(\xi)$. Critical points and values are illustrated in Appendix A. The following lemma expresses integral transforms as sums over critical points. It is key to our algorithms.

**Lemma 2.** Let $\varphi : \mathcal{P} \to \mathbb{Z}$ and let $\xi \in \mathbb{R}^n$ be generic on $\mathcal{P}$. Then,

$$ \mathcal{R}_{\varphi}(\xi, \cdot) = \sum_{x \in \mathrm{Crit}_{\varphi}^{-}(\xi)} \Delta_{\varphi}^{-}(\xi, x) \mathbf{1}_{\{\langle \xi; x \rangle\}} + \sum_{x \in \mathrm{Crit}_{\varphi}^{\mathrm{ord}}(\xi)} \Delta_{\varphi}^{\mathrm{ord}}(\xi, x) \mathbf{1}_{(\langle \xi; x \rangle, +\infty)}, $$

$$ \mathrm{ECT}_{\varphi}(\xi, \cdot) = \sum_{x \in \mathrm{Crit}_{\varphi}^{-}(\xi)} \Delta_{\varphi}^{-}(\xi, x) \mathbf{1}_{[\langle \xi; x \rangle, +\infty)}. $$

Furthermore, if $\kappa : \mathbb{R} \to \mathbb{C}$ is locally integrable and $\bar{\kappa}$ is a primitive of $\kappa$ on $\mathbb{R}$, one has:

$$ \mathrm{HT}_{\varphi}(\xi) = - \sum_{x \in \mathrm{Crit}_{\varphi}^{\mathrm{ord}}(\xi)} \Delta_{\varphi}^{\mathrm{ord}}(\xi, x) \cdot \bar{\kappa}(\langle \xi; x \rangle). $$
22:6 Efficient Computation of Topological Integral Transforms

Let $\mathcal{C}$ be a cubical complex. The result below can be adapted to arbitrary polytopal complexes but we do not make use of such a general statement in our algorithms. Since our definition of cubes only allow for axis-aligned cubes, any $\xi \in \mathbb{R}^n$ with non-zero coordinates is generic on $\mathcal{C}$. For any $\xi = (\xi_1, \dots, \xi_n) \in \mathbb{R}^n$ with non-zero coordinates we define its sign vector $\epsilon^\xi \in \{\pm 1\}^n$ by $\epsilon_i^\xi = \operatorname{sgn}(\xi_i)$. The following lemma ensures that only $2^n$ sets of critical points and values are necessary to compute integral transforms in practice.

**Lemma 3.** Let $\xi, \xi' \in (\mathbb{R} \setminus \{0\})^n$. If $\epsilon^\xi = \epsilon^{\xi'}$, then for any $\varphi : \mathcal{C} \to \mathbb{Z}$ and any $x \in \operatorname{Vert}(\mathcal{C})$ the sets and numbers $\operatorname{LwSt}(\xi, x)$, $\operatorname{UpSt}(\xi, x)$, $\Delta_{\varphi}^-(\xi, x)$, $\Delta_{\varphi}^{\operatorname{ord}}(\xi, x)$, $\operatorname{Crit}_{\varphi}^-(\xi)$ and $\operatorname{Crit}_{\varphi}^{\operatorname{ord}}(\xi)$ are respectively equal to the ones associated to $\xi'$.

**4 Implementation and asymptotic analysis**

The mathematical results of the previous section suggest algorithms allowing for efficient computation of the three integral transforms under consideration in this paper. In this section, we describe our implementation and provide worst-case time complexity bounds.

**4.1 Implementation of cubical complexes**

Our implementation of cubical complexes consists in a C++ class **EmbeddedComplex** that inherits from the cubical complex class **Bitmap_cubical_complex** from the GUDHI library (release 3.8.0) and additionally stores the embedding coordinates of the vertices in $\mathbb{R}^n$. This inheritance allows us to take any $n$-dimensional array/image as input and use GUDHI's initialization to turn it into a cubical complex. The GUDHI library stores the weighted cubical complex as a vector of bits representing the value assigned to each cell, with a specific convention for the indexing of cells. By default, one top-dimensional cell is associated to each pixel and the values of lower-dimensional cells are equal to the minimum of the values of their cofaces. By convention, the embedding coordinates of the vertices in $\mathbb{R}^n$ are evenly spaced and normalized in $[-0.5, 0.5]^n$. Storing the values of the complex in GUDHI's initialization is linear in the number of cells in the complex. We do not consider it as part of our algorithm. However, we take it into account in the running times in our experiments of Section 5.

**4.2 Preprocessing of critical points and values**

After the initialization, the first step to compute integral transforms using Lemma 2 is to compute the critical values $\Delta_{\varphi}^-(\epsilon, x)$ and $\Delta_{\varphi}^{\operatorname{ord}}(\epsilon, x)$ for all $x \in \operatorname{Vert}(\mathcal{C})$ and all sign vectors $\epsilon \in \{\pm 1\}^n$. While doing so, we store on the fly each one of the families $\Delta_{\varphi}^-(\epsilon) := \{\Delta_{\varphi}^-(\epsilon, x)\}_{x \in \operatorname{Crit}_{\varphi}^-(\epsilon)}$ and $\Delta_{\varphi}^{\operatorname{ord}}(\epsilon) := \{\Delta_{\varphi}^{\operatorname{ord}}(\epsilon, x)\}_{x \in \operatorname{Crit}_{\varphi}^{\operatorname{ord}}(\epsilon)}$ as a pair of two arrays, one for the critical points and one for the associated critical values, which induces no overhead in terms of asymptotic running time. Using GUDHI's implementation of weighted cubical complexes, accessing the dimension $\operatorname{dim}(\mathcal{C})$ and the value $\varphi(\mathcal{C})$ of a cell $\mathcal{C}$ takes constant time. Moreover, cells are indexed by a single number in this implementation. As a consequence, iterating through the $O(2^n)$ adjacent cells of a vertex requires arithmetic operations to compute indices of adjacent cells, resulting in $O(n2^n)$ time. The computation of critical values at a vertex thus has the same complexity, and the overall time complexity of computing critical values for all vertices and all sign vectors is then $O(n4^n \cdot \#\operatorname{Vert}(\mathcal{C}))$.
## 4.3 Exact representation and evaluation

After the preprocessing step, the computation of the transforms using Lemma 2 goes as follows. Let $\xi \in \mathbb{R}^n$, and denote $\#\text{Crit}_\varphi^-(\varepsilon^\xi)$ and $\#\text{Crit}_\varphi^{\text{ord}}(\varepsilon^\xi)$ respectively by $N^-$ and $N^{\text{ord}}$. First, we compute the sign vector $\varepsilon^\xi$ in $\mathcal{O}(n)$ time, which will always be negligible.

To compute the hybrid transform $\text{HT}_\varphi(\xi)$, we sum, over all ordinary critical points $x \in \text{Crit}_\varphi^{\text{ord}}(\varepsilon^\xi)$, the evaluation assumed to be done in constant time of the kernel $\bar{\kappa}$ on the dot products $\langle\xi; x\rangle$. Hence a running time in $\mathcal{O}(nN^{\text{ord}})$.

To compute the Euler characteristic transform $\text{ECT}_\varphi(\xi, \cdot)$, we first compute $T = \{\langle\xi; x\rangle\}_{x \in \text{Crit}_\varphi^-(\varepsilon^\xi)}$ in $\mathcal{O}(nN^-)$ time. Then, we sort the set $T$ and the sets $\text{Crit}_\varphi^-(\varepsilon^\xi)$ and $\Delta_\varphi^-(\varepsilon^\xi)$ by ascending values of $\langle\xi; x\rangle$ in $\mathcal{O}(N^- \log N^-)$ time. Finally, we compute a vector $E$ that contains the cumulative sums of the sorted values $\Delta_\varphi^-(\varepsilon^\xi)$ in $\mathcal{O}(N^-)$ time. An exact representation of the function $\text{ECT}_\varphi(\xi, \cdot)$ is then returned as the pair of vectors $(T, E)$. The total running time is in $\mathcal{O}(N^-(n + \log N^-))$.

To compute the Radon transform, we do the same as for the ECT but replacing $\text{Crit}_\varphi^-(\varepsilon^\xi)$ by $\text{Crit}_\varphi^{\text{ord}}(\varepsilon^\xi)$ and $\Delta_\varphi^-(\varepsilon^\xi)$ by $\Delta_\varphi^{\text{ord}}(\varepsilon^\xi)$, to get an analogously constructed pair $(T, E')$ in $\mathcal{O}(N^{\text{ord}}(n + \log N^{\text{ord}}))$ time. Then, we sort the elements $x \in \text{Crit}_\varphi^-(\varepsilon^\xi)$ by ascending values of $\langle\xi; x\rangle$ in $\mathcal{O}(N^-(n + \log N^-))$ time. After this, we use binary search to find, for each $x \in \text{Crit}_\varphi^-(\varepsilon^\xi)$, the critical point $y_x \in \text{Crit}_\varphi^{\text{ord}}(\varepsilon^\xi)$ with highest value $\langle\xi; y_x\rangle < \langle\xi; x\rangle$. This takes $\mathcal{O}(N^- \log N^{\text{ord}})$ time. We store $E'[y_x] + \Delta_\varphi^-(\varepsilon^\xi, x)$ for each $x \in \text{Crit}_\varphi^-(\varepsilon^\xi)$ in an array $E_c$. If there exist several points $x' \in \text{Crit}_\varphi^-(\xi)$ with the same dot product $\langle\xi; x'\rangle$, we only store the cumulative sum of the critical values $\Delta_\varphi^-(\varepsilon^\xi, x')$ in one entry of $E_c$. An exact representation of the function $\mathcal{R}_\varphi(\xi, \cdot)$ is then returned as the tuple of the four sorted vectors $T, E', T_c = \{\langle\xi; x\rangle\}_{x \in \text{Crit}_\varphi^-(\varepsilon^\xi)}$ and $E_c$, in $\mathcal{O}(N^{\text{ord}} \cdot n + (N^{\text{ord}} + N^-) \cdot \log N^{\text{ord}})$ time.

**Remark 4.** To better understand these complexities, note that the number of vertices of $\mathcal{C}$ is asymptotically equivalent to the number of pixels in the image. However, the total number of vertices appears only in the complexity of the preprocessing step. In many applications, the numbers $N^{\text{ord}}$ and $N^-$ are small compared to $\#\text{Vert}(\mathcal{C})$ so that the computation step is significantly optimized. This heuristic will be confirmed in Section 5.

To evaluate the Euler characteristic transform $\text{ECT}_\varphi(\xi, \cdot)$ at a point $t \in \mathbb{R}$ using its exact representation $(T, E)$, we use a binary search to find the element $x \in \text{Crit}_\varphi^{\text{ord}}(\varepsilon)$ with greatest $\langle\xi; x\rangle \le t$ and return $E[x]$ in $\mathcal{O}(\log N^{\text{ord}})$. We also implemented a vectorization routine to get the values of the ECT at sorted and evenly spaced points inside a specified interval. This is done by iterating simultaneously over $T$ and over the points of the interval, hence a running time in $\mathcal{O}(N + N^{\text{ord}})$.

The algorithm to evaluate the Radon transform using its exact representation $(T, E', T_c, E_c)$ is identical, except when there exists $x \in \text{Crit}_\varphi^-(\varepsilon)$ such that $t = \langle\xi; x\rangle$, in which case it returns $E_c[x]$. Again, we use binary search to determine the existence of such a point $x$ in $\mathcal{O}(\log N^-)$ time.

# 5 Experimental evaluation

In this section, we illustrate our optimizations on several batches of experiments. In all experiments, integral transforms are computed in 100 directions drawn uniformly at random in the cube $[0, 1]^n$. We compute hybrid transforms with $\kappa(t) = \sin(t)$. To allow for a comparison with demeter, we instantiate our cubical complexes using the dual construction, that is, initializing the values of the vertices with the pixel values in the images, and setting the values of higher-dimensional cells as the minimum of the values of their vertices. Moreover,
22:8 Efficient Computation of Topological Integral Transforms

we compare our exact version of the ECT to the discretized version of demeter with resolution T=100. This is justified by the fact that our vectorization step turning our exact version to a discretized one has negligible cost compared to the other steps, as shown in Figure 4.

Our code has been compiled with gcc 9.4.0 and run on a workstation with an Intel(R) Core(TM) i7-4770 CPU (8 cores, 3.4GHz) and 8 GB RAM, running GNU/Linux (Ubuntu 20.04.1). We measure time consumption using the Python module time. The input data and the benchmark Python scripts are available upon request.

5.1 Data set example

To begin with this experimental section, we compute integral transforms on a real-world data set: fashion_MNIST. This data set contains 60000 grayscale images of size 28 × 28 representing several types of clothes. We chose this data set because it contained global geometric information, namely the shape of clothes.

Results are shown in Table 1. In the first three lines, we timed our implementation of integral transforms on the fashion_MNIST data set with unchanged pixel values, potentially ranging from 0 to 255. In the next four lines, we timed our implementation as well as the ECT of demeter on a binarized version of the fashion_MNIST data set, where the 50% brightest pixels of each image are set to 1 and all other pixels to 0. This binarization was necessary since demeter does not cope with grayscale images.

<table><caption>Table 1 Timing (s) of our implementation and of the demeter ECT on the fashion_MNIST data set.</caption><thead><tr><th>Software</th><th>Init.</th><th>Preproc.</th><th>Comp.</th><th>Total</th></tr></thead><tbody><tr><td>Radon</td><td>13</td><td>143</td><td>290</td><td>447</td></tr><tr><td>HT</td><td>13</td><td>91</td><td>32</td><td>137</td></tr><tr><td>ECT</td><td>13</td><td>54</td><td>117</td><td>184</td></tr><tr><td>Radon (binary)</td><td>13</td><td>138</td><td>56</td><td>207</td></tr><tr><td>HT (binary)</td><td>13</td><td>88</td><td>14</td><td>115</td></tr><tr><td>ECT (binary)</td><td>13</td><td>52</td><td>32</td><td>97</td></tr><tr><td>demeter ECT</td><td>0</td><td>266</td><td>1544</td><td>1810</td></tr></tbody></table>

What stands out from this experiment is that all transforms can be computed in a reasonable amount of time on this large data set. The preprocessing step seems to be even more efficient in the case of hybrid transforms, where the evaluation step in 100 directions is done in between a third and a sixth of the time of the preprocessing step.

More importantly, this experiment perfectly illustrates the efficiency of our implementation. Each 28 × 28 image of the fashion_MNIST data set induces a cubical complex with 784 vertices. On average, there are 172 classical critical points with a standard deviation of 50 and 237 ordinary critical points with a standard deviation of 65. Thus, around a quarter of the vertices are classical critical points, and the same is true for ordinary critical points. In contrast, the complexes of the binarized data set have 23 classical critical points on average with a standard deviation of 12 and 36 ordinary critical points with a standard deviation of 19. As a consequence, after the initialization and preprocessing steps, our optimization using critical points allows us to compute an exact version of the ECT on the binarized version more than 50 times faster than demeter computes its discretized version. Overall, our implementation is 18 times faster than demeter.
<figure><img src="image_3.png" /><figcaption>(a) Initialization.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Preprocessing or “complexification” of demeter.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Computation times excluding initialization and preprocessing.</figcaption></figure>

Figure 1 Step-by-step computation times for varying size *m* of *m* × *m* images.

Moreover, thanks to our generalization of Euler critical values to weighted polytopal complexes, it is still 10 times faster to compute our exact version of the ECT on the original grayscale data set than to compute it on the binarized data set using demeter.

## 5.2 Asymptotic analysis

In this section, we perform an experimental asymptotic analysis of our algorithms. This analysis confirms our theoretical results on time complexities of Section 4.

### Size of images

We study the computation times on *m* × *m* images with respect to the parameter *m*. We use binary images with constant number of critical points made of a periodic pattern of black squares in a white background. The number of squares is constant when *m* varies to ensure a constant number of critical points. Namely, they are 100 classical critical points and 200 ordinary critical points for all images.

As one could expect, the initialization (Figure 1a) seems linear in the number of cells in the complex, that is, quadratic in *m*. Similarly, the preprocessing step (Figure 1b) is linear in the number of vertices in the complex, that is, quadratic in *m*. The strength of our implementation is illustrated by Figure 1c. After the initialization and the preprocessing steps, computing the integral transforms is done in constant time when the number of critical points is fixed. The curves for our times in Figure 1c are constant and lie between 1 and 5 milliseconds.

Overall, our computation times are dominated by the preprocessing step. Remarkably, our implementation has improved computation times by a factor of 4 compared with demeter.

### Number of critical points

We study the computation times with respect to the number of critical points in 1000 × 1000 images. We use the same periodic pattern of squares to control the number of critical points, this time with a varying number of squares.

The initialization and preprocessing steps are done in constant time for a constant size of image. More precisely the initialization takes around 1 second and the preprocessing of classical and ordinary critical points take approximately 3 and 6 seconds respectively. The “complexifying” step of demeter also takes around 6 seconds.
<figure><img src="image_6.png" /><figcaption>■ **Figure 2** Computation times for varying number of critical points excluding (2a) and including (2b) initialization and preprocessing steps.</figcaption></figure>

After the initialization and the preprocessing steps, computing the integral transforms is done in linear (HT) and linearithmic (ECT, Radon) times in the number of critical points for a fixed size of images. This is clearly illustrated in Figure 2a. Overall, computing an Euler characteristic curve for a binary image with four millions pixels takes between a 0.1 milliseconds to a few seconds depending on the number of critical points.

Interestingly, the computation time for the `demeter` package decreases with the number of critical points. This is due to our periodic patterns with varying size used to generate images with a varying number of critical points. The smaller the patterns are, the greater the number of critical points, but the lesser edges and two-dimensional cubes in the complex. As a consequence, when the number of critical points increases up to one quarter of the total number of vertices, `demeter` has a faster evaluation step than our evaluation of the ECT (Figure 2a). Time for the total pipeline are thus very similar between their and our implementation (Figure 2b).

However, this fact can be mitigated by two remarks. First of all, we compute an exact version of the Euler characteristic transforms and not an approximated version with a resolution of $T=100$. This approximated version is far from the exact one, as they are 25000 changes of values in the true Euler characteristic curve when they are 25000 classical critical points. Secondly, the number of critical points in real-world binary images containing global geometric features are likely to be way lower than a quarter of the vertices, as suggested by our study on a real-world data set in Section 5.1. While the number of critical points is likely to be higher for grayscale images, the `demeter` package can no longer be used in this case.

## Number of directions

We study the influence of the number of requested directions on the total running time. All computations are done on a random binarized image of the fashion MNIST dataset. We run the computation of the ECT (including initialization and preprocessing) for 1 to 10000 directions sampled uniformly in $[0, 1]^2$, and measure the corresponding times. The initial value of each curve corresponds to the initialization time. We give these plots on log-scale in Figure 3a, then on linear scale with the initial values forced to zero (which corresponds to ignoring the initialization time) in Figure 3b. By a linear regression, we check that the points we describe sample a linear function and we determine that its slope is $2.8 \cdot 10^{-6}$ seconds by direction for our method, and $2.3 \cdot 10^{-4}$ seconds by direction for `demeter`.
<figure><img src="image_7.png" /><figcaption>■ **Figure 3** Computation times (s) for varying number of directions of our and demeter's implementation of the ECT on a single image of the fashion-MNIST dataset.</figcaption></figure>

## Vectorization

On Figure 4 we ran the vectorization routine on 1000x1000 image with $10^6$ critical points for a varying dimension of vectorization, that is, sampling the ECT on an increasing number of points $t \in \mathbb{R}$. The complexity of the vectorization routine is confirmed by Figure 4. More importantly, vectorization is negligible (a few milliseconds) compared to other steps of the computation of the transform (a few seconds each).

<figure><img src="image_8.png" /><figcaption>■ **Figure 4** Vectorization time (s) with respect to dimension for a 1000 × 1000 image with $10^6$ critical points.</figcaption></figure>

## Average number of critical points

In this section, we expose an experimental study of the average number of critical points in random settings for two-dimensional grayscale images as the number of shades of gray increases. Pixel values are drawn using a uniform distribution in a first setting and using a Gaussian random field in a second setting.

In both cases, the number of both types of critical points in grayscale images is approximately half the number of vertices as soon as the number of pixel values is above 10. In the binary case, the number of critical points of each type is between 10% and 20% of the total number of vertices. These results corroborate the computation times observed on the data sets presented above: the number of critical points is in practice lower than the number of vertices, and much lower than the number of cells in the complex.
While this experiment shows that the number of critical point is of the same order as the number of vertices, it is expected not to be the case on digital images with global geometric features, as suggested in Section 5.1.

<figure><img src="image_9.png" /><figcaption>■ **Figure 5** Number of critical points in two random settings. Red and blue symbols represent respectively classical and ordinary critical points. Dots and crosses respectively represent these numbers for the uniform setting and for the Gaussian random field setting.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>■ **Figure 6** Relative pre-processing time (s) given the number of cores used for several complex sizes.</figcaption></figure>

## 5.3 Parallelization

We parallelized the routines `preproc_ordinaryiconv` and `preproc_classicaliconv` as they are the most time-consuming functions. To do so, we parallelized the inner loop on Vert($C$), each core handling the computations on a subset of the vertices. The readings are concurrent as cores only need to access fixed cell values. Each core returns an array of critical points and an array of critical values and the main core concatenates the results. Parallelizing the inner loop instead of the main one is efficient because the number of sign vectors in dimension $n$ is $2^n$, and $n = 2$ or 3 in most practical cases. We also implemented a parallelization of the routine `compute ht` that takes as input a vector of directions and returns a vector containing the values of the hybrid transforms in these directions. In this case each core handles the computations on a subset of the directions and returns the result to the main core that concatenates them. On Figure 6 we study the time of the preprocessing of critical points
(both classical and ordinary) given the number of cores used, relatively to the single core setup. We do this for several complex sizes with cells values uniformly sampled in {0,1}. We see that the parallelization induces a really small overhead, even for small complexes. Moreover, recall that even with only one core, our optimizations allow our software to be 100 times faster than the demeter package on images of the fashion_MNIST data set.

# 6 Conclusion

We presented a novel implementation of topological integral transforms on weighted cubical complexes, based on existing results from piecewise linear Morse theory slightly adapted to weighted polytopal complexes. These results allow for the fast computation of integral transforms on weighted cubical complexes. As a consequence, our implementation can efficiently compute exact versions of the ECT and Radon transform, and we showed experimentally that it is 18 times faster than the demeter ECT on a real-world data set, and 50 times faster after the preprocessing step. We also added a vectorization routine for data analysis purposes.

A natural extension of our implementation would be to arbitrary polytopal complexes. For any such complex, there is a cellular decomposition of (the dual of) $\mathbb{R}^n$ such that Lemma 3 holds for any $\xi, \xi'$ lying in the same cell [7]. This cellular decomposition is induced by the arrangement of hyperplanes that are orthogonal to the edges of the complex. The main challenge is that the size of this arrangement for a general polytopal complex with $n$ vertices in $\mathbb{R}^d$ may be as large as $n^{2d}$, as opposed to only $2^d$ for cubical complexes.

# References

1. Erik J Amézquita. demeter. https://github.com/amezqui3/demeter, 2021.

2. Erik J Amézquita, Michelle Y Quigley, Tim Ophelders, Jacob B Landis, Daniel Koenig, Elizabeth Munch, and Daniel H Chitwood. Measuring hidden phenotype: Quantifying the shape of barley seeds using the euler characteristic transform. in *silico Plants*, 4(1), 2022.

3. Mladen Bestvina and Noel Brady. Morse theory and finiteness properties of groups. *Inventiones mathematicae*, 129(3):445–470, 1997.

4. Leo Betthauser. Topological-Reconstruction. https://github.com/lbetthauser/Topological-Reconstruction, 2018.

5. Lorin Crawford, Anthea Monod, Andrew X. Chen, Sayan Mukherjee, and Raúl Rabadán. Predicting clinical outcomes in glioblastoma: An application of topological and functional data analysis. *Journal of the American Statistical Association*, 115(531):1139–1150, 2020.

6. Justin Curry, Robert Ghrist, and Michael Robinson. Euler calculus with applications to signals and sensing. In *Proceedings of Symposia in Applied Mathematics*, volume 70, pages 75–146, 2012.

7. Justin Curry, Sayan Mukherjee, and Katharine Turner. How many directions determine a shape and other sufficiency results for two topological transforms. *Transactions of the American Mathematical Society, Series B*, 9(32):1006–1043, 2022.

8. Robert Ghrist, Rachel Levanger, and Huy Mai. Persistent homology and euler integral transforms. *Journal of Applied and Computational Topology*, 2:55–60, 2018.

9. Robert Ghrist and Michael Robinson. Euler-Bessel and Euler-Fourier transforms. *Inverse Problems*, 27(12), 2011.

10. Olympio Hacquard and Vadim Lebovici. Euler characteristic tools for topological data analysis. arXiv preprint:2303.14040, 2023.

11. Vadim Lebovici. Hybrid transforms of constructible functions. *Foundations of Computational Mathematics*, pages 1–47, 2022.

12. Hugo Passe. eucalc. https://github.com/HugoPasse/Eucalc, 2021.
22:14 Efficient Computation of Topological Integral Transforms

13 Pierre Schapira. Operations on constructible functions. *Journal of Pure and Applied Algebra*, 72(1):83-93, 1991.

14 Pierre Schapira. Tomography of constructible functions. In *International Symposium on Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes*, pages 427-435. Springer, 1995.

15 Wai Shing Tang, Gabriel Monteiro da Silva, Henry Kirveslahti, Erin Skeens, Bibo Feng, Timothy Sudijono, Kevin K. Yang, Sayan Mukherjee, Brenda Rubenstein, and Lorin Crawford. A topological data analytic approach for discovering biophysical signatures in protein dynamics. *PLOS Computational Biology*, 18(5):1-42, May 2022.

16 Katharine Turner, Sayan Mukherjee, and Doug M. Boyer. Persistent homology transform for modeling shapes and surfaces. *Information and Inference: A Journal of the IMA*, 3(4):310-344, 2014.

17 Oleg Yanovich Viro. Some integral calculus based on euler characteristic. In *Topology and geometry—Rohlin seminar*, pages 127-138. Springer, 1988.

18 Günter M Ziegler. *Lectures on Polytopes, Updated seventh printing of the first edition*, volume 152. Springer, 2006.

A Running example

In this section, we present an example of run of our algorithms on a small binary cubical
complex. The complex is illustrated in Figure 7. Cells with non-zero values include 9 vertices,
9 edges, and 2 squares.

<figure><img src="image_11.png" /><figcaption>■ **Figure 7** Example of cubical complex. Colored cells have value 1 and light gray cells have value 0. The blue cells are vertices, the green ones are edges and the gray ones are faces. For greater clarity we display the keys of the vertices.</figcaption></figure>

A.1 Preprocessing

In this section, we unfold the computation of the critical points. We only detail it for the
ordinary ones, the computation for the classical ones being similar. We detail the operations
for ε = (1, -1) and x = 32. First we compute e-, the euler characteristic of the lower star
of x in direction ε. The coordinates of x are c = (32 ÷ 7, 32 mod 7) = (4, 4), so we explore
all cells with coordinates (4 - δ₀ε₀, 4 - δ₁ε₁) where (δ₀, δ₁) ∈ {0, 1}². These coordinates
are (4, 4), (3, 4), (4, 5) and (3, 5). They correspond to the cells 32, 31, 39 and 38 respectively.
Therefore, one has e_ = -1 + 1 + 0 + 0 = 0. In the same way we can compute e+, the upper
<table><caption>Table 2 Critical points and values of the complex in Figure 7.</caption><thead><tr><th>ε</th><th>Crit<sub>ϕ</sub><sup>ord</sup>(ε)<br>Δ<sub>ϕ</sub><sup>ord</sup>(ε)</th><th>Crit<sub>ϕ</sub><sup>-</sup>(ε)<br>Δ<sub>ϕ</sub><sup>-</sup>(ε)</th></tr></thead><tbody><tr><td>(-1, -1)</td><td>[28, 16, 44, 34]<br>[1, 1, -1, -1]</td><td>[0, 30, 44, 34]<br>[1, -1, 1, 1]</td></tr><tr><td>(-1, 1)</td><td>[42, 18, 32, 34]<br>[1, -1, 1, -1]</td><td>[0, 18, 32, 34]<br>[1, 1, -1, 1]</td></tr><tr><td>(1, -1)</td><td>[42, 18, 32, 34]<br>[-1, 1, -1, 1]</td><td>[0, 42]<br>[1, 1]</td></tr><tr><td>(1, 1)</td><td>[28, 16, 44, 34]<br>[-1, -1, 1, 1]</td><td>[0, 28, 16, 30]<br>[1, 1, 1, -1]</td></tr></tbody></table>

star at $x$ in direction $\epsilon$. We visit cells 32, 33, 25 and 26 and output $e_+ = 1 + -1 + -1 + 0 = -1$. Thus $e_+ \neq e_-$ and the vertex 32 is an ordinary critical point of $\epsilon = (1, -1)$ with critical value $-1$. We list critical points and values for all $\epsilon \in \{-1, 1\}^2$ in Table 2.

## A.2 Computation

Now, we unfold the computation step for Radon transform, ECT, and hybrid transform. Throughout the section, we set $\xi = (2, 2)$ and denote its sign vector simply by $\epsilon = (1, 1)$. We first describe the computation of the ECT. We start by computing the dot product of $\xi$ with all the critical points in $C_{\epsilon}^{-} = \{0, 28, 16, 30\}$. These critical points have coordinates $[-0.5, -0.5]$, $[-0.5, 0.167]$, $[-0.167, -0.167]$ and $[-0.167, 0.167]$ in $\mathbb{R}^2$. Thus the dot products are $[-2, 0.67, -0.67, 0]$. Then, we sort the critical points by increasing dot product, obtaining the list sorted($C_{\epsilon}^{-}$) = [0, 28, 16, 30]. Finally, we sum the critical values $\Delta_{\epsilon}^{-}$ in the order given by the sorting of critical points. We begin with $T = [-\infty]$ and $E = [0]$. After visiting the first two vertices, 0 and 28, we obtain $T = [-\infty, -2, -0.67]$ and $E = [0, 1, 2]$. The next sorted vertex is 16, which has the same dot product as 28. In that case, we add the variation associated to 16 to the last element of $E$. We obtain $T = [-\infty, -2, -0.67]$ and $E = [0, 1, 3]$. Finally, after adding vertex 30 we obtain $T = [-\infty, -2, -0.67, 0]$ and $E = [0, 1, 3, 2]$. The computation of the Radon transform follows the same process at first but with $C_{\epsilon}^{\text{ord}}$ instead of $C_{\epsilon}^{-}$. In that case, $T = [-\infty, -0.67, 0.67, 1.33]$ and $E' = [0.0, -2.0, -1.0, 0.0]$. For each $x \in C_{\epsilon}^{-}$, we must then find the greatest $i_0$ such that $T[i_0] \le \langle \xi; x \rangle$, and add $\langle \xi; x \rangle$ to $T_c$ and $E'[i_0] + \Delta_{\epsilon}^{-}[x]$ to $E_c$. We thus get $T_c = [-2.0, -0.67, 0.0]$ and $E_c = [1.0, 0.0, -1.0]$. Finally, computing the hybrid transform is done by summing $-\bar{\kappa}(\langle \xi; x \rangle) \cdot \Delta_{\epsilon}^{\text{ord}}[x]$ for $x \in C_{\epsilon}^{\text{ord}}$. In our case, with $\bar{\kappa} = \exp$, we obtain $\text{HT}_{\varphi}(\xi) \simeq -4.715$.

## A.3 Evaluation

We now unfold the evaluation of the Radon transform and ECT. In this example, we want to evaluate $\text{ECT}_{\varphi}(\xi, t)$ for some $t \in \mathbb{R}$. For $\xi = (2, 2)$, we computed in the above section that $T = [-\infty, -2, -0.67, 0]$ and $E = [0, 1, 3, 2]$. By a binary search, we find the index $i$ of the largest value in $T$ that is smaller than $t$. Then, we let $\text{ECT}_{\varphi}(\xi, t) = E[i]$. For instance $t = -0.3$ gives $i = 2$ and $t = 1$ gives $i = 3$. Thus, we get $\text{ECT}_{\varphi}(\xi, -0.3) = 3$ and $\text{ECT}_{\varphi}(\xi, 1) = 2$. The evaluation of the Radon transform is similar. First, we search for indices $i$ such that $T_c[i] = t$. In that specific case, the method outputs $\mathcal{R}_{\varphi}(\xi, t) = E_c[i]$. Apart from this, the computations are identical.
## A.4 Vectorization

The vectorization routine of the ECT (similar to the vectorization of the Radon transform) takes as input two bounds $a$ and $b$, and a number $N$ of points. Then, it returns $\text{ECT}_{\varphi}(\xi, a + \frac{i(b-a)}{N-1})$ for all $i \in [0, \dots, N-1]$. For instance, if $a = -2.0$, $b = 2.0$, and $N = 5$, we first iterate over $i = 0, \dots, 4$. For $i = 0$, we have $t_0 = a = -2$, so $\text{ECT}_{\varphi}(\xi, t_0) = E[1] = 1$. For $i = 1$, we know that $t_1 = -1 > t_0$, so we can start to search at index 1 in $T$. Then, we find that $\text{ECT}_{\varphi}(\xi, t_1) = E[1] = 1$. Finally, the values of $\text{ECT}_{\varphi}(\xi, t_0)$ for $t_0 \in [-2, -1, 0, 1, 2]$ are $[1, 1, 2, 2, 2]$.

# B Proofs of mathematical results

In this section, we prove Lemmas 2 and 3 on which our algorithms are built. We first state and prove some necessary lemmas before giving the proofs of our main results. Throughout the section, let $\mathcal{P}$ be a polytopal complex in $\mathbb{R}^n$ and let $\xi: \mathbb{R}^n \to \mathbb{R}$ be a linear form that is generic on $\mathcal{P}$.

These results and notions introduced in Sections 2 and 3 are more naturally phrased and proven using Euler calculus [13, 17], that is, the integral calculus of constructible functions with respect to the Euler characteristic. For instance, our definition of Euler characteristic $\chi(\varphi)$ of a weighted cubical complex $\varphi: \mathcal{P} \to \mathbb{Z}$ is equal to the Euler integral [13, Def. 2.2] of the constructible function $\hat{\varphi}: \mathbb{R}^n \to \mathbb{Z}$ defined by:

$$ \hat{\varphi} = \sum_{P \in \mathcal{P}} \varphi(P) \mathbf{1}_{\text{relint}(P)}, $$

where $\text{relint}(Q)$ denotes the relative interior of a polytope $Q \subseteq \mathbb{R}^n$. Similarly, for any $\xi \in \mathbb{R}^n$, the function $\mathcal{R}_{\varphi}(\xi, \cdot): \mathbb{R} \to \mathbb{Z}$ is equal to the so-called pushforward $\xi_*\hat{\varphi}$ in the language of Euler calculus; see [13, Eq. (2.5)]. For the sake of simplicity, we avoided the general definition of constructible functions and of Euler integration in the core of the paper. Along the section, we introduce notions of Euler calculus only when necessary to keep the text accessible to the broadest audience and refer to [13] for more details.

For any $P \in \mathcal{P}$ we denote by $\mathbf{1}_P: \mathcal{P} \to \mathbb{Z}$ the weighted cubical complex with value 1 on $P$ and 0 on any over cells (even on proper faces of $P$). We call $\mathbf{1}_P$ the indicating function of $P$. Then, any weighted cubical complex $\varphi: \mathcal{P} \to \mathbb{Z}$ can be written as a sum $\varphi = \sum_{P \in \mathcal{P}} \varphi(P) \mathbf{1}_P$.

The following lemma gives an explicit expression of the Radon transform of indicating functions. Before stating it, let us make the following remark. The Euler characteristic, Euler characteristic transforms, Radon transforms and hybrid transforms of a function $\varphi: \mathcal{P} \to \mathbb{Z}$ and of its restriction to the complex $\mathcal{P}^{\varphi}$ induced by the collection $\{P \in \mathcal{P}: \varphi(P) \neq 0\}$ are equal.

**Lemma 5.** Let $P \in \mathcal{P}$ and assume that $P$ is not a vertex of $\mathcal{P}$. Then, one has $\mathcal{R}_{\mathbf{1}_P}(\xi, \cdot) = (-1)^{\dim(P)-1} \mathbf{1}_{(\min_P f, \max_P f)}$.

**Proof.** This fact is an easy consequence of the Euler calculus formulation of our notions detailed in the beginning of this section. We give an elementary proof for completeness. Let us denote $\varphi = \mathbf{1}_P$. Since $\xi$ is linear, the set $P \cap \{\xi = t\}$ is a polytope (possibly empty) for any $t \in \mathbb{R}$. Moreover, since $\xi$ is generic on $\mathcal{P}$ and $P$ is not a vertex, this polytope is either a vertex of $P$ or it has dimension $\dim(P) - 1$.

Let $t \in \mathbb{R}$. We distinguish three cases. First assume that $t \notin [\min_P \xi, \max_P \xi]$. Since $P$ is connected and $\xi$ is continuous, one has that $[\min_P \xi, \max_P \xi] = \xi(P)$. Therefore, the polytope $P \cap \{\xi = t\}$ is empty, so is $\mathcal{P}^{\varphi}_{\{\xi=t\}}$ and hence $\mathcal{R}_{\varphi}(\xi, t) = 0$.
Assume now that $t \in \{\min_P \xi, \max_P \xi\}$. Since $\xi$ is linear and generic on $\mathcal{P}$, it must be that $P \cap \{\xi = t\}$ is a vertex $y$ of $P$. In that case $\mathcal{P}_{\{\xi=t\}}^{\varphi} = \{y\}$ and $\mathcal{R}_{\varphi}(\xi, t) = \varphi(y) = 0$ since $\varphi = \mathbf{1}_P$ and $P$ is not a vertex.

Suppose now that $t \in (\min_P \xi, \max_P \xi)$. In that case, the polytope $Q = P \cap \{\xi = t\}$ is neither empty nor a vertex. Therefore, this polytope has dimension $\dim(P) - 1$. Since $\mathcal{P}^{\varphi}$ consists of $P$ and all its faces, the polytopal complex $\mathcal{P}_{\{\xi=t\}}^{\varphi}$ contains only one polytope of dimension $\dim(P) - 1$ which is $Q$. Any other polytope $Q'$ of $\mathcal{P}_{\{\xi=t\}}^{\varphi}$ has dimension smaller or equal to $\dim(P) - 2$. In such a case, one has $\varphi_{\{\xi=t\}}(Q') \neq \varphi(P)$ and thus $\varphi_{\{\xi=t\}}(Q') = 0$. As a consequence, one has $\mathcal{R}_{\varphi}(\xi, t) = \varphi_{\{\xi=t\}}(Q)(-1)^{\dim(Q)} = (-1)^{\dim(P)-1}$. Hence the result. ◀

For any vertex $x$ of $\mathcal{P}$, we denote by $\mathrm{LwSt}^*(\xi, x)$ the set $\mathrm{LwSt}(\xi, x) \setminus \{x\}$. Similarly, we define the set $\mathrm{UpSt}^*(\xi, x)$. Moreover, for any function $h: \mathbb{R} \to \mathbb{R}$ we denote one-sided limits at any $t \in \mathbb{R}$ by:

$$h(t^{-}) = \lim_{s \to t} h(s) \quad \text{and} \quad h(t^{+}) = \lim_{s \to t} h(s),$$

whenever they are defined. Moreover, following the notations of the beginning of this section, Euler calculus ensures that for any $\varphi: \mathcal{P} \to \mathbb{Z}$ the function $\mathcal{R}_{\varphi}(\xi, \cdot) = \xi_* \hat{\varphi}$ is a finite sum of indicating functions of real intervals [13, Thm. 2.3.(i)], hence admits all one-sided limits both from above and from below. The following lemma provides additional details on the Radon transform of indicating functions. It will be crucial in our study of critical points (Lemma 7).

**Lemma 6.** Let $x \in \mathrm{Vert}(\mathcal{P})$ be such that $\langle \xi; x \rangle = t$ and let $P \in \mathrm{St}(x)$. One has

$$\mathcal{R}_{1_P}(\xi, t) = \begin{cases} 1 & (P = \{x\}), \\ 0 & (P \in \mathrm{LwSt}^*(\xi, x) \cup \mathrm{UpSt}^*(\xi, x)), \\ (-1)^{\dim(P)-1} & (P \notin \mathrm{LwSt}(\xi, x) \cup \mathrm{UpSt}(\xi, x)). \end{cases}$$

Moreover,

$$\mathcal{R}_{1_P}(\xi, t^{-}) = \begin{cases} 0 & (P \in \mathrm{UpSt}(\xi, x)), \\ (-1)^{\dim(P)-1} & (P \notin \mathrm{UpSt}(\xi, x)). \end{cases}$$

Similarly,

$$\mathcal{R}_{1_P}(\xi, t^{+}) = \begin{cases} 0 & (P \in \mathrm{LwSt}(\xi, x)), \\ (-1)^{\dim(P)-1} & (P \notin \mathrm{LwSt}(\xi, x)). \end{cases}$$

**Proof.** Throughout the proof, denote $\varphi = \mathbf{1}_P$. Let us first prove the first equality. The case $P = \{x\}$ is clear. If $P \in \mathrm{LwSt}^*(x) \cup \mathrm{UpSt}^*(x)$, then $\mathcal{P}_{\{\xi=t\}}^{\varphi} = \{x\}$ and $\mathcal{R}_{\varphi}(\xi, t) = \varphi(x) = 0$. Now, suppose that $P \notin \mathrm{LwSt}(\xi, x) \cup \mathrm{UpSt}(\xi, x)$. Since $t \in \xi(P)$, the polytope $Q = P \cap \{\xi = t\}$ is non-empty. Moreover, the polytope $Q$ is not a vertex, for otherwise the generic linear form $\xi$ would be constant on the edge between the distinct vertices $Q$ and $x$. Therefore, the polytope $Q$ has dimension $\dim(P) - 1$. As in the proof Lemma 5, one can thus show that $\mathcal{R}_{\varphi}(\xi, t) = \varphi_{\{\xi=t\}}(Q)(-1)^{\dim(Q)} = (-1)^{\dim(P)-1}$.

Now, let us prove the second equality. If $P \in \mathrm{UpSt}(\xi, x)$, then $P \subseteq \{\xi \ge t\}$ and $\{\xi = s\} \cap P = \emptyset$ for all $s < t$. Hence $\mathcal{R}_{\varphi}(\xi, t^{-}) = 0$. If $P \in \mathrm{St}(x) \setminus \mathrm{UpSt}(\xi, x)$, then $\min_P \xi < t$. Since both $\min_P \xi$ and $t$ lie in the image of the connected set $P$ by the continuous map $\xi$, then $[\min_P \xi, t] \subseteq \xi(P)$. Therefore, for any $s \in [\min_P \xi, t]$, the polytope
22:18 Efficient Computation of Topological Integral Transforms

$Q_s = \{\xi = s\} \cap P$ is non-empty. The same argument as in the previous paragraph shows that $Q_s$ is not a vertex and hence has dimension $\dim(P) - 1$. Once again, we can conclude that $\mathcal{R}_{\varphi}(\xi, t) = \varphi_{\{\xi=t\}}(Q_s)(-1)^{\dim(Q_s)} = (-1)^{\dim(P)-1}$.

The proof of the last equality follows from the last one applied to the linear form $-\xi$. ◀

**► Lemma 7.** Let $\varphi : \mathcal{P} \to \mathbb{Z}$. For any $t \in \mathbb{R}$,

$$ \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \sum_{\substack{x \in \text{Crit}_{\varphi}^{\neg}(\xi) \\ \langle \xi; x \rangle = t}} \Delta_{\varphi}^{\neg}(\xi, x), $$

$$ \mathcal{R}_{\varphi}(\xi, t^{+}) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \sum_{\substack{x \in \text{Crit}_{\varphi}^{\text{ord}}(\xi) \\ \langle \xi; x \rangle = t}} \Delta_{\varphi}^{\text{ord}}(\xi, x). $$

**Proof.** Let us prove the first equality. Let $t \in \mathbb{R}$. Since $\varphi = \sum_{P \in \mathcal{P}} \varphi(P) \mathbf{1}_P$ one has:

$$ \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \sum_{P \in \mathcal{P}} \varphi(P) (\mathcal{R}_{\mathbf{1}_P}(\xi, t) - \mathcal{R}_{\mathbf{1}_P}(\xi, t^{-})). \quad (B.1) $$

Fix $P \in \mathcal{P}$. Suppose that no vertex of $P$ belongs to the level set $\{\xi = t\}$. Since $\xi$ is linear, it implies that $t \neq \min_P f$ and $t \neq \max_P f$. Thus $t \in (\min_P f, \max_P f)$ or $t \in \mathbb{R} \setminus [\min_P f, \max_P f]$ and Lemma 5 implies that:

$$ \begin{aligned} \mathcal{R}_{\mathbf{1}_P}(\xi, t) - \mathcal{R}_{\mathbf{1}_P}(\xi, t^{-}) &= (-1)^{\dim(P)-1} \left( \mathbf{1}_{(\min_P \xi, \max_P \xi)}(t) - \mathbf{1}_{(\min_P \xi, \max_P \xi)}(t^{-}) \right) \\ &= 0. \end{aligned} $$

Now, suppose there exists a vertex $x$ of $P$ such that $\langle \xi; x \rangle = t$. Lemma 6 implies that if $P \not\in \text{LwSt}(\xi, x)$, then $\mathcal{R}_{\mathbf{1}_P}(\xi, t) - \mathcal{R}_{\mathbf{1}_P}(\xi, t^{-}) = 0$. Indeed, either $P \in \text{UpSt}^*(\xi, x)$ and then $\mathcal{R}_{\mathbf{1}_P}(\xi, t) = \mathcal{R}_{\mathbf{1}_P}(\xi, t^{-}) = 0$, or $P \not\in \text{LwSt}(\xi, x) \cup \text{UpSt}(\xi, x)$ and then $\mathcal{R}_{\mathbf{1}_P}(\xi, t) = \mathcal{R}_{\mathbf{1}_P}(\xi, t^{-}) = (-1)^{\dim(P)-1}$. Therefore, (B.1) can be written:

$$ \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \sum_{\substack{x \in \text{Vert}(P) \\ \langle \xi; x \rangle = t}} \sum_{P \in \text{LwSt}(\xi, x)} \varphi(P) \left( \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) \right). \quad (B.2) $$

where the sum in the right-hand side is over all polytopes $P \in \mathcal{P}$ such that there exists $x \in \text{Vert}(P)$ with $\{\xi = t\}$ and $P \in \text{LwSt}(\xi, x)$. Moreover, for any $P \in \mathcal{P}$ there is at most one vertex $x$ of $P$ such that $P \in \text{LwSt}(\xi, x)$. Indeed, if there were $x \neq x'$ such that $P \in \text{LwSt}(\xi, x) \cap \text{LwSt}(\xi, x')$, then $\langle \xi; x \rangle = \max_P \xi = \langle \xi; x' \rangle$. Denote $m = \langle \xi; x \rangle = \langle \xi; x' \rangle$. Then $\{\xi = m\} \cap P$ is a face of $P$ that is not a vertex and on which the generic linear form $\xi$ is constant, a contradiction.

Therefore, (B.2) can be written:

$$ \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \sum_{\substack{x \in \text{Vert}(P) \\ \langle \xi; x \rangle = t}} \sum_{P \in \text{LwSt}(\xi, x)} \varphi(P) \left( \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) \right). \quad (B.3) $$

Let $x \in \text{Vert}(P)$ such that $\langle \xi; x \rangle = t$ and $P \in \text{LwSt}(\xi, x)$. Lemma 6 implies that $\mathcal{R}_{\mathbf{1}_P}(\xi, t) = 1$ if $P = \{x\}$ and $0$ if $P \in \text{LwSt}^*(\xi, x)$. The same lemma implies that $\mathcal{R}_{\mathbf{1}_P}(\xi, t^{-}) = 0$ if $P = \{x\}$ and $(-1)^{\dim(P)-1}$ if $P \in \text{LwSt}^*(\xi, x)$. Therefore, one has:

$$ \begin{aligned} \sum_{P \in \text{LwSt}(\xi, x)} \varphi(P) \left( \mathcal{R}_{\mathbf{1}_P}(\xi, t) - \mathcal{R}_{\mathbf{1}_P}(\xi, t^{-}) \right) &= \varphi(x) - \sum_{P \in \text{LwSt}^*(\xi, x)} (-1)^{\dim(P)-1} \varphi(P) \\ &= \sum_{P \in \text{LwSt}(\xi, x)} (-1)^{\dim(P)} \varphi(P) \\ &= \Delta_{\varphi}^{\neg}(\xi, x). \end{aligned} $$
Therefore, (B.3) can be written:

$$ \mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \sum_{\substack{x \in \text{Vert}(\mathcal{P}) \\ \langle \xi; x \rangle = t}} \Delta_{\varphi}^{-}(\xi, x). \qquad (B.4) $$

The set $\text{Crit}_{\varphi}^{-}(\xi)$ being the set of vertices $x$ with non-zero $\Delta_{\varphi}^{-}(\xi, x)$, the first equality of the lemma follows.

Now, let us prove the second equality. Replacing $\xi$ by $-\xi$, one gets from (B.4) that:

$$ \mathcal{R}_{\varphi}(\xi, t^{+}) - \mathcal{R}_{\varphi}(\xi, t) = \sum_{\substack{x \in \text{Vert}(\mathcal{P}) \\ \langle \xi; x \rangle = t}} \chi(\varphi|_{\text{UpSt}(\xi, x)}). \qquad (B.5) $$

Hence the result, using (B.4), (B.5), and the fact that $\mathcal{R}_{\varphi}(\xi, t^{+}) - \mathcal{R}_{\varphi}(\xi, t^{-}) = \mathcal{R}_{\varphi}(\xi, t^{+}) - \mathcal{R}_{\varphi}(\xi, t) - (\mathcal{R}_{\varphi}(\xi, t) - \mathcal{R}_{\varphi}(\xi, t^{-}))$. $\blacksquare$

## B.1 Proof of Lemma 2

Let us prove the equality for the Radon transform. Euler calculus [13, Thm. 2.3.(i)] ensures that $\xi_{*}\tilde{\varphi}$ is a finite sum of indicating functions of real intervals. Therefore, it is locally constant on the complement of a finite number of points in $\mathbb{R}$. It is then an easy exercise to check that Lemma 7 yields the result:

$$ \mathcal{R}_{\varphi}(\xi, \cdot) = \sum_{x \in \text{Crit}_{\varphi}^{-}(\xi)} \Delta_{\varphi}^{-}(\xi, x) \cdot \mathbf{1}_{\{(\xi; x)\}} + \sum_{x \in \text{Crit}_{\varphi}^{\text{ord}}(\xi)} \Delta_{\varphi}^{\text{ord}}(\xi, x) \cdot \mathbf{1}_{((\xi; x), +\infty)}. \qquad (B.6) $$

By definition of the hybrid transform, integrating (B.6) against a locally integrable kernel $\kappa: \mathbb{R} \to \mathbb{C}$ yields the result for $\text{HT}_{\varphi}(\xi)$.

The formula for the ECT can be easily shown using the classical convolution operation of Euler calculus [13, Sec. 4]. This bilinear operator takes two functions $\theta$ and $\theta'$ which are finite sums of indicating functions of intervals of $\mathbb{R}$ and returns a function $\theta \star \theta'$ with the same property. One can show that $\text{ECT}_{\varphi}(\xi, t) = \mathcal{R}_{\varphi}(\xi, \cdot) \star \mathbf{1}_{[0, +\infty)}(t)$; see for instance [11, Ex. 7.4, Prop. 7.5]. The expression for the ECT follows then from (B.6) and the fact that for $a \in \mathbb{R}$ one has $\mathbf{1}_{(a, +\infty)} \star \mathbf{1}_{[0, +\infty)} = 0$ and $\mathbf{1}_{a} \star \mathbf{1}_{[0, +\infty)} = \mathbf{1}_{[a, +\infty)}$; see [13, Exs. 4.1, 4.3]. This concludes the proof of Lemma 2. $\blacksquare$

## B.2 Proof of Lemma 3

Let $\mathcal{C}$ be a cubical complex. In that case, any edge $\{x, y\} \in \mathcal{C}$ is such that $x - y$ is colinear to a canonical basis vector of $\mathbb{R}^n$. Let $\xi$ and $\xi'$ be in $\mathbb{R}^n \setminus \{0\}$ with $\varepsilon\xi = \varepsilon\xi'$. Since $\xi$ and $\xi'$ have same sign vector, we have $\text{sgn}\langle\xi; x - y\rangle = \text{sgn}\langle\xi'; x - y\rangle$ for all edge $\{x, y\} \in \mathcal{C}$. It is an easy exercise to check that thus $\xi$ and $\xi'$ induce the same order on $\text{Vert}(\mathcal{C})$. As a consequence, the lower and upper stars of vertices of $\mathcal{C}$ with respect to $\xi$ and $\xi'$ are the same. Moreover, the sets and numbers $\Delta_{\varphi}^{-}(\xi, x)$, $\Delta_{\varphi}^{\text{ord}}(\xi, x)$, $\text{Crit}_{\varphi}^{-}(\xi)$ and $\text{Crit}_{\varphi}^{\text{ord}}(\xi)$ depend only on $\text{LwSt}(\xi, x)$ and $\text{UpSt}(\xi, x)$ (and not on $\xi$), hence the result. $\blacksquare$