DIGITAL EULER CHARACTERISTIC TRANSFORM

HENRY KIRVESLAHTI AND XIAOHAN WANG

arXiv:2411.08522v1 [math.AT] 13 Nov 2024

Abstract. The Euler Characteristic Transform (ECT) of Turner et al. pro-

vides a way to statistically analyze non-diffeomorphic shapes without relying
on landmarks. In applications, this transform is typically approximated by a
discrete set of directions and heights, which results in potential loss of infor-
mation as well as problems in inverting the transform. In this work we present
a fully digital algorithm for computing the ECT exactly, up to computer pre-
cision; we introduce the Ectoplasm package that implements this algorithm,
and we demonstrate this is fast and convenient enough to compute distances in
real life data sets. We also discuss the implications of this algorithm to related
problems in shape analysis, such as shape inversion and sub-shape selection.
We also show a proof-of-concept application for solving the shape alignment
problem with gradient descent and adaptive grid search, which are two pow-
erful methods neither of which is possible using the discretized transform.

1. Introduction
Shapes are difficult to analyze statistically, because they have no obvious metric
structure needed for statistical machinery. Without a metric structure, we can not
define such basic statistical quantities like mean of two shapes. A standard idea
in shape analysis is to pick points on the shapes, declare correspondence between
them, and define the distance between shapes to be the distance between the ma-
trices defined by the landmarks, modulo appropriate Euclidean actions. This leads
to what is called the Kendall Shape Space. Using landmarks leads to loss of in-
formation, because shapes are not point clouds, and there is no canonical choice
for picking the correspondence between the landmarks. This is especially relevant
for the modern age, where our shape data are represented in high fidelity in online
repositories.
Another idea is that of diffeomorphism-based shape alignment. This corresponds
to a digitalization of Kendall Shape Space as in this approach the shapes are rep-
resented as surfaces and their distances are studied via correspondence maps, con-
tinuous objects that can achieve arbitrary precision. A key restriction for there
to be such a correspondence map is that the shapes must be diffeomorphic. This
assumption is not realistic for a multitude of shape data, such as skulls that may
have different number of cavities or medical data such as brain MRI.
A third idea for shape analysis is based on integral geometry and geometric
measure theory. A powerful tool in this is the Euler Characteristic Transform
(ECT), which is mathematically a topological Radon transform. This transform

2020 Mathematics Subject Classification. Primary: 62R40, 68U05; Secondary: 55N31.

Key words and phrases. Topological data analysis, Statistical Shape Analysis, Topological
Radon Transform.
1
2 HENRY KIRVESLAHTI AND XIAOHAN WANG

is the subject of this work. This framework allows for statistical analysis beyond
landmarks without relying on the diffeomorphism assumption.
All existing methods for computing the ECT rely on evaluating the ECT on a
discrete grid of directions, and most of the time, also the heights. This is somewhat
undesirable if our goal is to define a digital tool for shape analysis - such an approach
still requires us to pick a discrete set of directions, and our inference is based on
these, a setup much similar as with landmarks. While there is a theoretical notion
of a Moduli space for shapes that can be told apart from each other by finitely
many directions [5], this number is very large, indeed much larger than what is
used in practice. Also, the result states that such transform is injective for these
shapes, but there is to date even a theoretical guarantee on how many directions
give a good notion of distance for a given collection of shapes. In practice, this
can be only solved by picking more and more directions until the distances seem to
stabilize.
In this work we introduce the Ectoplasm, a fully-digital shape analysis tool
for dissimilar shapes. The work relies heavily on the observations made in [5] and
folklore results regarding the analytic form of the transform. We push these results
further to develop a fully digital ECT algorithm for shape analysis. In some regard,
this work complements the theoretical approach to unify shape theory via sheafifi-
cation presented in [1], with the distinction that from a theoretical perspective the
Persistent Homology Transform (PHT) has nicer properties than the ECT, which
is preferable from applied perspective and subject of this work.
The main advantages of this framework are threefold: First, the downstream
analysis, for example statistical inference, is improved, because the digital trans-
form makes use of all the data. This is especially relevant for getting accurate
information on which shapes are closest to each other, a key ingredient in mani-
fold learning applications that rely on small distances to approximate the tangent
space. Second, not having to discretize the transform obviates the problem of hav-
ing to choose the discretization parameters, a major headache for the practitioners.
Third, as the transform is computed in its exact form, the rich mathematical struc-
ture of the transform can be exploited to its full extent. For example, the digital
algorithm allows for a genuine O(d)-invariant action. We will discuss these dis-
cretization related issues and ways to circumvent them in Section 2. Finally, the
digital approach also gives us a ground truth, which can be used to assess how
accurate a given discretization is.
The paper is organized as follows: in Section 2 we review some of the background
on ECT and its typical use cases. In Section 3 we give a simplified, brute-force ver-
sion of the digitalization algorithm in two dimensions.This serves as an introduction
to Section 4, where we present the 3D-algorithm. The technical details pertaining
to spherical integrals are in A. In Section 5 we apply the the algorithm to perform
exploratory data analysis on a real life data set, and compare the results against the
conventional, discretization based method, as well as previously published methods
based on landmarks and diffeomorphisms. In Section 6, we study the two aforemen-
tioned solutions to the alignment problem that the digital transformation allows
for.
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 3

2. Background on the ECT

In this section we provide a brief background on the ECT and its use cases
in applications. We will take as concrete approach as possible when introducing
the background material, a more mathematically oriented reader can consult, for
example [5] for a more general treatment of the subject. All our objects are finite
and compactly supported.
The applied upshot of the ECT is that it allows us to compare shapes. This
statement necessitates the question of what constitutes a shape. In general, the
ECT framework is defined in the flexible language of o-minimal structures, but as
this work is mainly concerned with, and indeed only applicable to, piecewise linear
meshes, we will state the relevant results in this language. This is also the relevant
framework for practical considerations, because computer standard architecture
prefers triangles over curved geometric objects.
Our meshes will be finite geometric simplicial complexes, with two meshes de-
clared identical if they are identical as subsets of Rd . We will recall the somewhat
technical definition of a geometric simplicial complex.
Definition
  2.1. A (finite) abstract simplicial complex Σ is a union of subsets of
N (called simplices) that is closed under taking subsets, i.e. if τ ≤ σ and σ ∈ Σ,
then τ ∈ Σ. The singletons of Σ are called its vertices, and the dimension of a
simplex is its cardinality minus one.
Definition 2.2. Let Σ be a finite abstract simplicial complex with vertices x1 , . . . , xn
and M simplices. A geometric realization map of Σ is an affine map f from a subset
of Rn to Rd (for d some positive integer), whose domain consists of
(1) The coordinate vectors ei ;
(2) Any strictly convex combination of a collection of coordinate vectors
{ej1 , ej2 , . . . , ejm }, whenever {xj1 , xj2 , . . . , xjm } is a simplex of Σ.
Note that by affinity, the geometric realization is completely determined by the
image of the coordinate vectors ei , and the combinatorial structure of Σ.
Definition 2.3. A geometric simplicial complex X is an abstract simplicial complex
Σ, together with an injective geometric realization map f .
For our purposes, this definition is too fine, because we are only interested in
meshes as subsets of the Euclidean space, typically R3 . Hence we don’t care about
the combinatorial structure beyond that. This conundrum was also discussed in
[5], Remark 5.12.
Definition 2.4. A mesh is a geometric simplicial complex as a subset of Rd .
We will use the word mesh to distinguish it from a shape, which is usually re-
garded as a mesh, modulo certain Euclidean symmetries, such as rotations, trans-
lations and scalings.
We will require our meshes to be supported in the unit ball. This restriction,
albeit unsatisfactory from the mathematical perspective, is well-grounded in the
context of Kendall shape space and the requirement that scalar multiples of shapes
are equal, as well as the practical wisdom that a collection of meshes should be
scaled to the same measure before they are analyzed.
One of the important invariants of a simplicial complex X is its Euler Charac-
teristic. It is indeed a topological invariant of the complex, but for our purposes
we can take it to be the alternating sum of its simplex counts.
4 HENRY KIRVESLAHTI AND XIAOHAN WANG

Definition 2.5. The Euler Characteristic χ(X) of a geometric simplicial com-

plex X with b0 0-simplices (vertices), b1 1-simplices (edges), and in general bj
j-simplices, j = 0, . . . , n, is
n
X
χ(X) = (−1)j bj .
j=0

The following fact is a consequence of the o-minimal Triangulation Theorem

from [13].
Proposition 2.6. If mesh M is obtained from two different geometric simplicial
complexes X and Y , then χ(X) = χ(Y ), so that χ(M ) is well-defined.
Definition 2.7. For a mesh X ⊂ Rd , and a direction v ∈ S d−1 , we can subset X
at height h:

Xv (h) = {x ∈ X|x · v ≤ h}.

We will call the subset Xv (h) the filtration of X to height h at direction v.
For a fixed v, this is actually a filtration of topological spaces indexed by h ∈ R
in that Xv (h1 ) ⊆ Xv (h2 ) whenever h1 ≤ h2 .
Note that the restriction {x ∈ Rd |x · v ≤ h, x ∈ X}, is a piecewise linear set,
but it need not be a subcomplex of X. However, we can effectively treat it as such,
which is the subject of the following definition and theorem.
Definition 2.8. Write Xv,h for the subcomplex of X consisting of the simplices
that are contained in the half-space
Hv,h = {x ∈ Rd |x · v ≤ h}.
We call Xv,h the restriction of X to height h at direction v.
We are now ready to define the Euler Characteristic Transform, the main subject
of this work:
Definition 2.9. Let X ⊂ Rd be a mesh, and Xv,h the subcomplex from Definition
2.8. The Euler Characteristic Transform ECTX of X is an integer-valued function
on S d−1 × R defined as

ECTX (v, h) = χ(Xv,h ).

While we don’t need the general definition for Euler Characteristic of a topo-
logical space for computing anything, we note the following result which is a con-
sequence of [2], and it justifies the filtration viewpoint of Definition 2.9. We state
this result for reference:
Theorem 2.10. For a mesh X ⊂ Rd , and a direction v ∈ S d−1 , and height h ∈ R,
we have:
χ(Xv (h)) = χ(Xv,h ).
The upshot of a more theoretical framework is that by a well-known consequence
of the work by Schapira [11] is that ECTX (v, h) determines X completely. This is
also relevant for what is done in this paper, as the results that we need follow from
this framework. However, we do not need this machinery for the exposition of the
work that we do here, readers are referred to [5] for the details.
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 5

2.1. The ECT distance. With the ECT, we can define the ECT inner product
between two meshes X and Y as

Z 1 Z 
(1) ⟨X, Y ⟩ = ECTX (v, h)ECTY (v, h) dvdh .
−1 S d−1
This of course gives us a notion of norm, or distance between two meshes:
(2) ||X − Y ||2 = ⟨X − Y, X − Y ⟩ = ⟨X, X⟩ − 2⟨X, Y ⟩ + ⟨Y, Y ⟩.
The expression (2) of squared distance above is especially convenient for our
purposes. We note of the equivalent expression, which is more common in the
literature:

Z 1 Z

2 1/2
(3) dECT (X, Y ) = ECTX (v, h) − ECTY (v, h) dvdh .
−1 S d−1

Computing the distances in Equation (3) seems like a daunting task, for the
integral is taken over the sphere. In practice, this expression is always approximated
via some sort of discretization. While the integral indeed is very complicated, with
the help of a computer it can be computed analytically. An effective algorithm for
doing so is the main result of this paper. Before introducing this algorithm, we will
make a few comments on the discretization and its properties.
2.2. Discretization of ECT. A straightforward way to discretize the transform
and to make Equation (3) computable is to pick a set of direction v1 , v2 , . . . , vn
and a set of heights h1 , . . . , hm and then represent the (discretized) transform as
an n × m vector (or matrix). For these discretized objects, (3) is nothing more
than the L2 norm of the matrices. If we hook this idea back to landmarks, we are
ultimately performing a similar analysis: Our inference is based on matrices, that
represent finitely many directions evaluated at finitely many heights. This bears
resemblance to the landmark problem - how do we choose the directions? We will
discuss the merits and demerits of the approximation in this subsection.
One of the main results in this direction is the complexity class result in [5],
Theorem 7.14. – that there is an upper bound on how many directions are required
to tell two shapes in a certain complexity class apart. This is a mathematical
result, meaning that by picking that many directions guarantees that the discretized
transform is injective on the complexity class. This result, being one of the only
ones in the literature, is important, but we don’t have sufficient understanding on
what is the best procedure to pick the discretization in practice for a collection of
shapes. We can not, for example, make statements that N many directions explain
99% of the variability in distance, or to say which directions should be chosen for
best representation of the shapes if we can only afford n many directions. The lack
of these type of results makes using these methods harder, because there is no clear
cut way to choose the discretization parameters.
This problem is exacerbated when the statistical question is something more
complicated than regression or classification. If our only goal is to do classifica-
tion, and we use a kernel method, we only need to worry about approximating
the distance matrix well, because all of our inference is based on that. Then the
conventional wisdom to pick n and m as large as possible, and pick the directions
(and heights) as uniformly as possible, seems to be a reasonable advice. On the
6 HENRY KIRVESLAHTI AND XIAOHAN WANG

other hand, if we are interested in performing a more complicated task, such as fea-
ture selection, which was done in [14], then the choice of discretization parameters
becomes a double-edged sword. One one hand, we would like to increase n and m
to make the distances more accurate. But on the other hand, if we pick n and m
too large, the variation in the n × m regression parameters can too homogeneous,
causing at best multicollinearity, and at worst, making the parameters not even
estimable.
The digitalization procedure proposed in this work alleviates these two prob-
lems: By not using a discretization, the distances are computed (up to computer
precision) in closed form, so we know our inference is as good as it can get. Also,
the digitalization obviates the need to choose the discretization parameters in the
first place, so that multicollinearity of features will not be an issue. Of course, one
needs to rework all the statistical machinery for variable selection, because with
the digital transform we replace vectors with functions, so that the statistical ma-
chinery must involve functional data analysis. This interesting problem is left for
future research.
Another advantage that the digital version has is that it is better suited for
capturing small differences in the meshes. This is especially relevant for applica-
tions that make use of tangent space approximations, such as ISOMAP. Because
smaller distances are used for constructing the tangent space, a method that is more
accurate in getting the small distances right will give more accurate inference.
Furthermore, the discretized version suffers a small inconvenience when it comes
equivariance. Mathematically, the ECT is an O(d)-equivariant transform, meaning
that for a rotation matrix (possibly including reflection) A, we have

ECTAX (v, h) = ECTX (Av, h).

If one uses the discretized version, one only has access to finitely many direc-
tions. In a typical use case, these directions do not form an orbit under A, so that
we don’t have access ECTX (Av, h) from ECTX (v, h). This means they must be
somehow approximated, and this complicates the analysis and can lead to costly
computations. See [14] Section 4 in the supplementary material for details.
On the other hand, from the digitalized version, ECTX (Av, h) is readily available
from ECTX (v, h). What is more, the alignment penalty map

A 7→ −2⟨X, AY ⟩
is, for generic meshes X and Y , almost everywhere differentiable. This phenom-
enon has been discussed in e.g. [9]. This means that in theory, with the digital
transform, the alignment problem is amenable to a gradient descent type solution
(rather than grid search). We study this further in Section 6.
Finally, the discretized version has no hope of using the Schapira’s inversion
formula to go from transforms back to shapes. With the digitalized version of the
transform, it is conceivable that Schapira-type inversion formula could be used for
pulling back evidence back on the original problem.

3. The Algorithm
The purpose of this Section is to explain how we can evaluate Expression (3) ana-
lytically. To do so, we first introduce a brute force algorithm for the 2D case, where
the formulae are simpler and the transform is easier to visualize. We postpone the
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 7

actual 3D-algorithm with additional refinements to Section 4. The algorithm, while

quite complicated, makes use of several well-known facts of the ECT. We will state
these facts shortly. With these facts, we will describe the transforms (or equiv-
alently the meshes) in what we will call the proto-transfom format. This format
is a bridge between the mesh and its transform, and a key to the digitalization
algorithm.
The first well-known fact about the ECT that we want to state is that for a
piecewise linear mesh, the Euler Curve fv (h) can only change at the vertices.
Proposition 3.1. Let X be a piecewise linear mesh with vertices {x1 , . . . , xn } and
v a direction such that ⟨xi , v⟩ =
̸ ⟨xj , v⟩ for i ̸= j and fv (h) the Euler curve in
direction v. Then fv (h) can only change at the vertices, that is at heights hi given
by hi = ⟨v, xi ⟩ for some vertex xi of X.
Definition 3.2. Following the notation and assumptions of Proposition 3.1 , vertex
xi is called Euler-critical (with respect to v) and hi is called a critical height. We
will call

ai (v) := fv (hi ) − fv (hi − ϵ),

where ϵ is a small positive number, for example minj,k ⟨xj − xk , v⟩/2,
the gain associated to vertex xi in direction v.
Next we describe a consequence of the stratification result in [5]. The argument
stated here has also appeared previously in [12].
Proposition 3.3. Let X be a mesh with vertices {x1 , . . . , xn }. If two directions v
and ν satisfy

⟨v, xi ⟩ ≤ ⟨v, xj ⟩ if and only if ⟨ν, xi ⟩ ≤ ⟨ν, xj ⟩

for all i, j ∈ [n], the directions v and ν are equivalent, and they have the same
critical vertices with same gains.
Proofs of Propositions 3.1, 3.3. See Propositions 5.18 and Lemma 5.19 of [5]. □
Next we will state a simple fact of the height function with an explicit formula
in dimension 2.
Proposition 3.4. Let p = (x, y) be a vertex
in R2 . Its height function h with
respect to v parametrized as v = cos(t), sin(t) is given by
hp (v) = p · v = x cos(t) + y sin(t),
so that its integral over a polygon Pk = [τ, θ] is given by:

Z θ
IPk (p) = cos(t)x + sin(t)y dt
τ   
= cos(τ ) − cos(θ) y + sin(θ) − sin(τ ) x.

Proof. A straightforward calculation. □

With these facts we can now define the proto-transform format, as well as a
brute force algorithm for computing it in two dimensions. The idea behind the
8 HENRY KIRVESLAHTI AND XIAOHAN WANG

proto-transform format is to list the finite pieces of the transform in a systematic

format that allows us to compute Expression (3) with ease.
Let X be a 2D-mesh with vertices {x1 , . . . , xn }.
• Solve ⟨v, xi ⟩ = ⟨v, xj ⟩ for all i ̸= j, get vi,j , −vi,j
·
• Order the vi,j by their angle, get pieces Pk = [θk , θk+1 ], k = 1, . . . , N .
θk +θk+1
• Set αk = 2
• Look at the order of xi under ⟨αk , xi ⟩.
• For each xi , evaluate the discrete ECT at height hi + ϵ, i.e. just above xi
(with e.g. ϵ from Definition 3.2). Record this value to ak,i , and record xi
to bk,i .
SN
• The ECT of X is now represented as k=1 (Pk , ak , bk ).
SN
We call the the collection k=1 (Pk , ak , bk ) the proto-transform format of X.
Next we will show how we can use this format to compute integrals of the form 3
in 2-dimensions.
The first observation is that inside the pieces Pk , the order of the vertices stays
the same. Then, inside each Pk (and heights [−1, 1]), the ECTX (v, h) is a simple
function described as
nk
X

ECTX|Pk (v, h) = ai,k 1≥⟨v,bk,i ⟩ (h)1≤⟨v,bk,i+1 ⟩ (h) ,
i=1
where we define ⟨v, bk,nk +1 ⟩ to be identically 1.
This expression is again nothing new, it is equivalent to the formula in Lemma
5.19 of [5], expressed in explicit coordinates.
This is a simple function of the height functions from Proposition 3.4, whose
integral is readily available.
To compute the (squared) distance in Equation (2), we first note that each of
the three pieces

dECT (X, Y )2 = ||X − Y ||2 = ⟨X, X⟩ − 2⟨X, Y ⟩ + ⟨Y, Y ⟩.

are just simple functions of the height functions from Proposition 3.4, because
the product of simple functions is simple.
SN
Explicitly, if the ECT of X is represented as k=1 (Pk , ak , bk ) and we write
Pk = [θk , θk+1 ], the first term ⟨X, X⟩ can be expressed as:

Z 1 Z X nk
M X
a2i,k 1≥⟨v,bk,i ⟩ (h)1≤⟨v,bk,i+1 ⟩ (h) dv dh


⟨X, X⟩ =
−1 S d−1 k=1 i=1
M Z
X 1 Z nk
X
a2i,k 1≥⟨v,bk,i ⟩ (h)1≤⟨v,bk,i+1 ⟩ (h) dv dh



=
k=1 −1 Pk i=1
M  nk
X X

= a2nk ,k (θk+1 − θk ) − a21,k IPk (b1,k ) + IPk (bi,k ) a2i+1,k − a2i,k .
k=1 i=2

Here IPk (bi,k ) is the integral over the height function of the vertex bi,k over
polygon Pk , which was explicitly computed in Proposition 3.4. The integral is then
just a sum of trigonometric functions, and analytically tractable.
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 9

Similar reasoning applies for the other two terms, but note that the cross-product
term in the middle requires us to overlay the pieces Pk of X and Qk Y to a common
subdivision. See Example 3.5 for an explicit computation.
Note that the algorithm that we have described in this section is just for illustra-
tion purposes. In fact, it is very inefficient. For

two meshes with n and m vertices,
we would first need to compute n2 and m


2 pairwise intersections to get Pk , Qk ,
and then we still have to perform the expensive operation of overlaying these two
triangulations. Luckily, this complexity can be greatly reduced. As a matter of fact,
we don’t have to compute all the pairwise intersections, and we can also make the
pieces Pk larger than what was described here. We will discuss these improvements
in the next section, where we also introduce the algorithm for 3 dimensions, which
is also more interesting from the perspective of real-life applications.

Vertex heights against angle

4
2
height

0
−2
−4

0 1 2 3 4 5 6

angle

Figure 1. Example 3.5, the height functions of the 4 vertices. The

vertical lines represent equi-height angles, or endpoints of each Pk .
Vertices 1,2,3 and 4 are represented in black, red, green and blue,
respectively.

Example 3.5. In this example we compute the ECT inner product of a 2D-triangle
by hand, following the brute force algorithm. Let Σ be an abstract simplicial complex
{(123, 234)}, and consider its geometric realization X given by map fX : e1 7→
(−2, −1), e2 7→ (0, 1), e3 7→ (0, 4) and e4 7→ (2, 0). For simplicity of computations,
we will not scale X to the unit ball (nor center it), instead we shall carry the
integration out over the origin-centered ball of radius 4.
We first proceed by solving the equi-height angles. For two pairs
of vertices
−x2
(x1 , y1 ) and (x2 , y2 ), these are given by θ = atan(D) = atan xy21 −y 1
These are
given in the Table 1 below, and also depicted in Figure 1.
Ordering the angles from table 1, and investigating the order bk of the vertices
xi inside each Pk together with the Euler characteristic ak , we can represent the
ECT of X in a table (Table 2), with Pk = [θk , θk+1 ]:
10 HENRY KIRVESLAHTI AND XIAOHAN WANG

Table 1. Auxiliary computations for solving equi-height angles

x1 y1 x2 y2 D Name
−2 −1 0 1 −1 D12
−2 −1 0 4 − 25 D13
−2 −1 2 0 −4 D14
0 1 0 4 0 D23
0 1 2 0 2 D24
1
0 4 2 0 2 D34

Table 2. The proto-transform format of mesh X. Each row lists

a polygon, together with the order of vertices inside the polygon
(bk ), and (ak ) lists the Euler characteristic after each vertex. For
convenience, the angles are expressed in terms of trigonometric
functions.

k tan(θk ) tan(θk+1 ) cos(θk ) sin(θk ) bk ak

1 0 1/2 1√ 0√ (p1 , p2 , p3 , p4 ) (1, 1, 1, 1)
2 1/2 2 2/√5 1/√5 (p1 , p2 , p4 , p3 ) (1, 1, 1, 1)
3 2 −4 1/ √5 2/√ 5 (p1 , p4 , p2 , p3 ) (1, 2, 1, 1)
4 −4 −1 −1/ √17 4/ √17 (p4 , p1 , p2 , p3 ) (1, 2, 1, 1)
5 −1 −2/5 −1/√ 2 1/√ 2 (p4 , p2 , p1 , p3 ) (1, 1, 1, 1)
6 −2/5 0 −5/ 29 2/ 29 (p4 , p2 , p3 , p1 ) (1, 1, 1, 1)
7 0 1/2 −1√ 0√ (p4 , p3 , p2 , p1 ) (1, 1, 1, 1)
8 1/2 2 −2/√5 −1/√5 (p3 , p4 , p2 , p1 ) (1, 1, 1, 1)
9 2 −4 −1/√ 5 −2/√ 5 (p3 , p2 , p4 , p1 ) (1, 1, 1, 1)
10 −4 −1 1/ √17 −4/ √17 (p3 , p2 , p1 , p4 ) (1, 1, 1, 1)
11 −1 −2/5 1/√ 2 −1/√ 2 (p3 , p1 , p2 , p4 ) (1, 1, 1, 1)
12 −2/5 0 5/ 29 −2/ 29 (p1 , p3 , p2 , p4 ) (1, 1, 1, 1)

Table 2 contains all the information to compute ⟨X, X⟩. This decomposition
contains a lot of polygons and the expression could be much simplified. This ineffi-
ciency is a consequence of using the brute force algorithm. In fact, we see that the
integrand is 1 except for parts between p2 and p4 in polygon P3 and between p2 and
p1 in polygon P4 .
From Table 2, we can readily compute ⟨X, X⟩. The integrand is 1 from the
minimum vertex until height 4, except for the two polygons P3 and P4 where the
squared Euler characteristic is 4 after p1 (in P3 ) and p4 (in P4 ) up until p2 . In
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 11

other words:
Z 4 Z X 4
12 X
a2i,k 1≥⟨v,bk,i ⟩ (h)1≤⟨v,bk,i+1 ⟩ (h) dv dh


⟨X, X⟩ =
−4 S d−1 k=1 i=1
12 Z θk+1 h Z θ3 Z θ4
X


i
= 4 − min(pi · v) dv + 4 p2 · v − p4 · v dv + p2 · v − p1 · v dv
θk i θ2 θ3
k=1
h i
= 4 2π + IP3 (p2 ) + IP4 (p2 ) − IP3 (p4 ) − IP4 (p1 )
X 7
X 11
X
− IPk (p1 ) − IPk (p4 ) − IPk (p3 )
k={1,2,3,12} k=4 k=8
1 1
1 1
4 2

= 8π + 4 √ + √ +4 − √ + √ −8 √ − √
5 17 17 2 17 5
1 4
−1 1
9 9 2 8 8 20
+8 √ − √ +4 √ + √ + √ + √ + √ + √ + √ + √
2 17 17 2 17 29 5 17 5 29
√
√ √
√
= 8 π + 2 + 3 2 5 − 17 + 29
≈ 42.878.

4. 3D Algorithm
In this Section, we describe an algorithm for computing Expression (2) for 3D-
meshes X and Y . The algorithm proceeds as follows: Input: Mesh X with vertices
{x1 , . . . , xn }.
(1) For each vertex xi , compute its star si .
(2) For each star si , compute the triple intersection points; that is, for all
xj ∈ si , xk ∈ si , xj ̸= xk ̸= xi ̸= xj , solve v · xi = v · xj = v · xk , get a
collection of antipodal points on S d−1 , one pair for each triplet. Label each
such point with the triplet (i, j, k) that generated it
(3) Triangulate S d−1 by placing the vertices with points computed from (2).
Connect two vertices if they share two of their labels and the inner products
of the points they define is non-negative.
(4) For each triangle in the triangulation, evaluate gain in the local ECT in the
star. If this is zero, discard the triangle, otherwise record the change.
(5) Merge adjacent triangles, that is, triangles that share an edge and that have
the same gain. This results in a spherical polygon.
(6) Get: For each vertex xi a collection of polygons Pk with their local gains
ak .
This results in having the mesh X in the proto-transform format, i.e. it is repre-
sented as an indicator function of height functions over polygons, which this time,
may be overlapping. After the algorithm, the next step of this algorithm is to use
this representation to compute the distances. This resembles by and large the case
introduced in the previous section, with relevant details about spherical integra-
tion left into the Appendix A. We will use the rest of section to give mathematical
reasoning behind Algorithm 2 as well as discuss some improvement ideas.
The first difference to the brute force algorithm is that we are now computing
the polygon intersections inside each link. When mesh resolution is increased, n,
the number of vertices, increases, but the star of each vertex will have more or less
the same number of vertices, for
example for primate molars k ≈
6. Using the links
we need to compute about n · k3 intersections, instead of the n3 intersections that
the brute force algorithm requires.
12 HENRY KIRVESLAHTI AND XIAOHAN WANG

This is based on the following, well-known, observation that allows one to com-
pute changes in Euler Characteristic solely based on local information:
Proposition 4.1. Let X be a piecewise linear mesh with vertices {x1 , . . . , xn } and
v a direction such that hi := ⟨xi , v⟩ =
̸ ⟨xj , v⟩ for i ̸= j and fv (h) the Euler curve
in direction v. Let also ϵ = mini̸=j |(hi − hj )|/2. The gain

fv (hi + ϵ) − fv (hi − ϵ)
associated to vertex xi in the Euler Curve fv is completely determined by Star(xi )
and v.
Proof. The Euler characteristic of X is the alternating sum of its simplex counts.
Any new simplices that are added to X when xi is added to it are by definition
contained in the star of xi . □
The second difference that the algorithm has is that now the spherical polygons
may or may not be disjoint. This is not much of a problem, because regardless of
which decomposition we use, we have the following:
Proposition 4.2.PnSuppose shapes K1 , K2 are
Pm compactly-supported on unit ball B d ,
let ECT(K1 ) = i=1 αi fi and ECT(K2 ) = i=1 βj gj , we denote αn+j = −βj and
fn+j = gj for simplicity, then we have
Z Z
2
d (K1 , K2 ) = (ECT(K1 ) − ECT(K2 ))2 dtdσ
S2 I
Z Z X n m
X
= ( αi fi − βj gj )2 dtdσ
S2 I i=1 j=1
Z Z m+n
X
= ( αi fi )2 dtdσ
S2 I i=1
m+n
X Z Z
= αi αj fi fj dtdσ
i,j=1 S2 I

m+n
X Z Z
= αi αj 1{Pi ∩Pj } · 1[max{x(Pi )·v, x(Pj )·v},1] dtdσ
i,j=1 S2 I

m+n
X Z
= αi αj 1 − max{x(Pi ) · v, x(Pj ) · v} dtdσ
i,j=1 Pi ∩Pj

m+n
X 
= αi αj Vol(Pi ∩ Pj )
i,j=1
Z
− x(Pi ) · v dtdσ
Pi ∩Pj ∩{x(Pi )·v>x(Pj )·v}
Z 
− x(Pj ) · v dtdσ .
Pi ∩Pj ∩{x(Pi )·v>x(Pj )·v}

To derive the closed formula of ECT distance, it suffices to derive the closed
formula of integral of x · v on arbitrary spherical polygons, and a formula for the
volume of the intersection of the polygons. In general, it is difficult to provide a
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 13

closed-form formula for the former spherical integral in d dimensions. The simple
case of 2 dimensions was already addressed in the previous section (Proposition
3.4), and the three-dimensional formula, which is the most important one for appli-
cations, is presented in Appendix A. Note that in three dimensions the intersection
area can be computed with, for example, the Gauss-Bonnet Theorem.

5. Case Study
Now that we have described our algorithm, we demonstrate its utility by applying
it to a practical problem. We also compare it to the distances obtained from
landmarks as well as a diffeomorphism based method [3]. The comparison may
be of independent interest to theoreticians pondering the differences between the
continuous Procrustes distances and the one obtained from ECT, and to what
extent these two are commensurate.

5.1. The Data. The dataset comprises 116 primate molars, representing diverse
genera and dietary preferences. This dataset was first analyzed in [3].

5.1.1. Data Processing. The meshes are aligned based on landmarks using the
Auto3DGM[4, 10] algorithm, with the furthest point sampling and 128 and 256
landmarks for phases 1 and 2, respectively.
For the purpose of the ECT methods, the meshes were converted to OFF files
with meshlab, and further subsampled to 1,000 faces with its quadric edge collapse
decimation algorithm. We further center and scale the meshes to the unit ball by
a standard procedure.

5.2. Methods. We compute the Digital ECT distances based on the algorithm de-
scribed in this manuscript, using the software provided in Section 7. The discretized
ECT distances are based on 326 directions obtained from a regular subdivision of
a regular octahedron. This is a standard procedure to produce close to uniformly
spread points on the 2-sphere. For each direction, the discrete Euler Characteristic
Transform is evaluated with 100 equally spaced heights between -1 and 1, resulting
in a matrix of size 362 × 100.
The landmark distance is based on the 2nd phase of the Auto3DGM algorithm
that was used to align the meshes. This algorithm makes use of minimum spanning
tree to propagate distances, unlike our ECT distances. The continuous procrustes
distance is a verbatim copy of the results in [3], reproduced from its accompanying
code, provided here for reference.

5.3. Results. To see a high-level overview of the different distances and how they
cluster the data, multidimensional Scaling Plots of the distances are presented in
Figure 2. These are 2-D summaries of the distance matrix, projecting points to
the plane such that the global distance is distorted as little as possible. From the
plots we see that all the methods give, on high level, similar kind of clustering of
the dietary groups.
In Table 3 The mantel correlation of the distance matrices are presented. These
are standard measures of correlations between distance matrices, accounting for all
the pairwise distances. We see that all the four methods are related, as the MDS
plots in Figure 2 hint. Unsurprisingly, we see that the Digital ECT is closely related
to the Discrete one, indicating that the discretization that we chose is quite good
14 HENRY KIRVESLAHTI AND XIAOHAN WANG

(a) Landmarks (b) Continuous Procrustes Distance

(c) Discrete ECT (d) Digital ECT

Figure 2. Multidimensional Scaling plots of the different distances.

The clustering is quite similar across the methods, with Omnivores and
Frugivores forming their own clusters.

at capturing the distances as a whole. Interestingly, the digital method has smaller
correlation with the two non-ECT methods than the discrete one.
In Figure 3 we also look at the t-sne plots of the collection computed with differ-
ent notions of distance. The t-sne, unlike MDS, puts emphasis on small distances,
which is a key element in, for example, manifold learning and Gaussian process
classifiers, that make use of local structure of the collection in their inference. It
is conceivable that there may be relevant small differences that yield to different
inference even if the distances are globally in great agreement. While we see some
differences between the two methods, the closest neighbors are in much agreement,
for example, the frugivores are far from insectivores, and the folivores form three
clear subclusters whose relative positions are similar in either display.
In spite of the small differences visible in t-sne plots, there is little difference in
the clustering behavior, and indeed if we tried to predict, for example, frugivores
against other classes using a local smoothing type method, we would get similar
inference.
The conclusion of this case study is that the discretization yields results that
are very accurate and the inference is for all intents and purposes just as good as
with the digital transform. This corroborates the rule of thumb that 326 uniformly
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 15

(a) Digital ECT (b) Discrete ECT

Figure 3. t-sne plots based on the distances from the two ECT meth-
ods.

Table 3. Correlation of pairwise distances (the Mantel correla-

tions) of the distance matrices produced by different methods.

Method 1 Method 2 Correlation

Digital ECT Discrete ECT 0.9991
Digital ECT Landmarks 0.6066
Digital ECT Continuous Procrustes 0.4896
Discrete ECT Landmarks 0.6070
Discrete ECT Continuous Procrustes 0.4927
Landmarks Continuous Procrustes 0.7373

spread directions with hundred sublevel sets is typically a good enough discretiza-
tion for a collection of surface-like 3D-shapes. In general, the procedure outlined
in this section can be used to objectively assess how good a given discretization is.
One such assessment is presented in Table 4, where we report the Mantel correla-
tions against the digital algorithm against various discretizations. Ultimately, the
discretization must be chosen in terms of whatever problem is studied, what we
present here is but a single indication of the discretization error on a high level.

Table 4. The Mantel correlations of the distances from different

discretization against the true values obtained from the digital al-
gorithm. The discretizations are obtained by subsetting the full
discretization of 100 sublevelsets and 326 directions uniformly.

Directions \ Sublevel sets 100 50 25 20 10 5

326 0.9991 0.9989 0.998 0.9973 0.9923 0.9740
110 0.9979 0.9975 0.9945 0.9923 0.9782 0.9460
38 0.9921 0.9912 0.9838 0.9786 0.9360 0.8481
14 0.9535 0.9487 0.9308 0.9142 0.8173 0.7077
6 0.8715 0.8610 0.8153 0.7921 0.6639 0.5413
16 HENRY KIRVESLAHTI AND XIAOHAN WANG

6. Shape Alignment with the Digital Transform

In general, given a collection

of shapes X1 , . . . , Xn , the shape alignment problem
refers to finding a set of n2 rigid transformations called correspondences ci,j that
best align each shape Xj to another shape Xi . This is quite a nuanced problem,
because the composition of two correspondences ci,j ◦ ck,i may not necessarily yield
the correspondence between Xj and Xk , so some sort of harmonization procedure
is in order, for example through a minimum spanning tree. In this section we study
the more modest problem of aligning two shapes X and Y , which is a first step to
studying the more general alignment problem. Solving the more general problem
efficiently is left for future.
Concretely, the problem that we want to solve is to find A ∈ SO(3) that mini-
mizes the distance

Z 1 Z

2 1/2
(4) dECT (X, AY ) = ECTX (v, h) − ECTAY (v, h) dvdh .
−1 S d−1

This is of course the same as maximizing the generalized cross-correlation

Z 1 Z
⟨X, AY ⟩ = ECTX (v, h)ECTAY (v, h) dvdh.
−1 S d−1

One of the immediate advantages that the digital transform has vis-à-vis the dis-
crete one is that it is automatically SO(d)-equivariant, and it is almost everywhere
differentiable.
If we are using the discrete transform, we have restricted the ECT to a finite
set of directions (v1 , . . . vn ). In particular, we do not have access to the ECT at
directions (Av1 , . . . Avn ), so we are required to approximate the action of A by a
permutation of the directions (v1 , . . . vn ). For example, in [14] this problem was
solved by discretizing SO(3) to a grid of size 16,000, which were approximated by
a permutations of the set of 2,918 directions. For each of the 16,000 actions Ai , we
would then solve the linear assignment problem

n
X
p̃i = arg minp:[n]→[n] (pvj − Ai vj )2 ,
j=1

which is a computationally very expensive operation.

In this Section we propose two algorithms for aligning shapes with Ectoplasm:
adaptive grid search and gradient descent. These two methods can be used in
tandem. These methods have three advantages over the existing ECT-based align-
ment schemes. Firstly, it removes the computational overhead of approximating
rigid motions with permutations. Secondly, the method involves no permutations,
removing an additional source of noise. Thirdly, we are not restricted to any par-
ticular a priori chosen grid, so we can simply add more iterations as we go to
gain further accuracy without having to redefine the grid and without solving the
expensive permutation approximation problem from scratch.
We will set Y be a mesh and X = AY with a randomly chosen A ∈ SO(3) so
that we have a well-defined minimum and we can study the convergence of our
algorithms.
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 17

6.1. Adaptive Grid Search. We begin by discretizing the group of rotations

SO(3). The group of rotations can be parametrized by 3 angles, (α, β, γ), which
can be seen as the quotient space of R3 with the following identifications:

(α, β, γ) ≃ (α + 2π, β, γ)
(α, β, γ) ≃ (α, β + 2π, γ)
(α, β, γ) ≃ (α+, β, γ + 2π)
(α, β, γ) ≃ (α + π, β + π, γ + π).
The angles α, β, γ are the axis-aligned rotations, and we can associate to each a
rotation matrix for example as follows:

   
cos α − sin α 0 cos β 0 sin β 1 0 0
Aα,β,γ =  sin α cos α 0  0 1 0  0 cos γ − sin γ  .
0 0 1 − sin β 0 cos β 0 sin γ cos γ
Next we pick an initial grid to start our search. We choose a grid of 8×8×5 = 320
points. Concretely, the grid is the cCrtesian product of the following values for
parameters α, β and γ:

α = 0, π/4, π/2, 3π/4, . . . , 7π/4;

β = 0, π/4, π/2, 3π/4, . . . , 7π/4;
γ = 0, π/4, π/2, 3π/4, π.
After each iteration, we pick the set of angles that produced the maximal cross-
correlation. For each element of this maximal angle, we take its two neighbors, and
subdivide the resulting interval to 4 pieces. Our new grid is the Cartesian product
of these intervals.
This way, each angle is on the same scale, and our grid is of size 43 = 64, a
multiple of 16, leading to efficient parallelization.
As a demonstration of this algorithm, we take as Y one of the shapes analyzed in
the previous section, and we take a random element A from SO(3) and set X = AY .
Using the notation above, in terms of (α, β, γ), A is represented as approximately
(4.76, 1.07, 2.10).
We present results of this run in Table 5. For each iteration i, we report the
angles, the ECT inner product, as well as the SO(3) distance between A and its
approximation Ai .
Recall the SO(3) distance can be computed for two rotation matrices A and Ai
as

trace(AATi ) − 1
d(A, Ai ) = arccos .
2
From Table 5 we see that the inner product is inversely proportional to the
SO(3) Distance, so that the inner product is in practice a useful metric for the
shape alignment problem, and the algorithm does not, for example, seek a local
minimum. The method seems to be able to find the correct alignments using
relatively fast.
18 HENRY KIRVESLAHTI AND XIAOHAN WANG

Table 5. Results of the Adaptive Grid Search

Iteration α β γ ⟨X, Ai Y ⟩ SO(3) Distance

1 4.71239 0.78540 1.57080 22.11458 0.56745
2 4.97419 1.04720 2.35619 23.54182 0.12827
3 4.79966 1.22173 2.18166 23.20256 0.15730
4 4.68330 1.10538 2.06531 24.01233 0.05861
5 4.76087 1.02781 2.14288 23.90973 0.06097
6 4.81258 1.07952 2.19459 24.14080 0.05507
7 4.77811 1.04504 2.09116 24.22715 0.03779
8 4.75512 1.06803 2.11415 24.39189 0.01881
9 4.73980 1.08335 2.06818 24.37854 0.02024
10 4.75002 1.07313 2.09882 24.49974 0.00863
11 4.75683 1.06632 2.09201 24.51984 0.00783
Truth 4.75875 1.07199 2.09902 24.60740 0.00000

6.2. Gradient Descent. The digital version of the shape alignment problem is, in
theory, almost differentiable, and the piecewise components of it are easily described
with even pen and paper. However, writing the gradient down exactly is quite a
complicated task. To this end, we also implement an auto-differentiable version
of the code, which allows us to compute the gradient by abusing the simplicity of
individual operators in conjunction with the chain rule. In order to do this, we
provide a proof-of-concept implementation of the algorithm with torch functions,
a popular machine-learning framework that has built-in automatic differentials for
basic operations. Due to the cost of the computational resources required for an
efficient implementation, we demonstrate a proof-of-concept of this procedure by
applying the algorithm on a simple non-convex polyhedron on 9 vertices and 14
faces obtained from simplifying a shape studied in the previous section. This is a
mesh with high curvature.
The ECT representation introduced in this paper and implemented in the code
package Ectoplasm allows for straightforward implementation of both vanilla gra-
dient descent and also stochastic gradient descent. Pn
In terms of this
Pn representation, write ECTX (v, h) = i=1 αi P (xi , vi ) and
ECTY (v, h) = j=1 αj P (yj , vj ), the alignment problem corresponds to finding the
optimal A that maximizes

Z 1 Z
⟨X, AY ⟩ = ECTX (v, h)ECTAY (v, h) dvdh
−1 S d−1
Z 1 Z n
X m
X
= αi P (xi , vi ) αj P (Ayj , Avj ) dvdh.
−1 S d−1 i=1 j=1

The gradient of this expression with respect to the vector angles (α, β, γ) is then
readily available through auto-differentiation, and we can seek the optimum via
gradient ascent.
For the stochastic gradient, we have an additional choice to make. We could
divide up the double sum over one index or the other, or even both of them. The
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 19

latter is not necessarily a good option1, because then we would at each iteration
compute the gradient based on the intersection of two sparse sets of spherical poly-
gons on S 2 . The gradient would then point out to whichever direction maximizes
the overlap between these two sets. This type of sampling would lead to large
volatility that would override the signal. By subsetting only one index, we ensure
the other set of polygons is full, i.e., there is a polygon (possibly several) in each
direction, so that the intersection cannot be maximized solely by seeking maximal
overlap.
We seek the alignment by vanilla gradient ascent, which we initialize with angles
θ0 := (α, β, γ) = (0, 0, 0). We update the angles θi of alignment Ai with the gradient
ascent

∇θi
θi+1 = θi + λ ,
||∇θi ||
with λ = 1 for the first 30 iterations and 0.1 after the 20 following iterations and
0.01 thereafter. The results are presented in Table 6.

Table 6. Results of the Gradient Ascent

Iteration α β γ ⟨X, Ai Y ⟩ SO(3) Distance

0 0 0 0 29.4329 2.90493
5 1.2458 0.3611 0.8870 29.8906 2.59675
10 0.6660 0.8395 0.4906 29.8233 2.84400
15 1.2837 1.0522 0.1899 29.7568 2.02619
20 2.6830 0.9540 -0.9083 31.1978 0.37889
25 1.8636 1.0757 -0.3600 30.4473 0.94800
30 2.6630 0.9739 -0.9524 31.1811 0.39757
35 2.3444 1.0060 -0.7960 32.0700 0.06131
40 2.4250 1.0194 -0.8537 32.0291 0.07331
45 2.3589 0.9934 -0.7833 32.0694 0.06037
50 2.4251 1.0178 -0.8535 32.0321 0.07292
55 2.3914 1.0061 -0.8184 32.2554 0.00587
60 2.3894 1.0097 -0.8083 32.2477 0.00717
65 2.3911 0.9995 -0.8162 32.2475 0.00726
70 2.3861 1.0089 -0.8149 32.2600 0.00398
Truth 2.3886 1.0059 -0.8151 32.2767 0.00000

While the gradient ascent seems to find the correct alignment without very com-
plicated optimization of λ, it does suffer from the fact that there are many local
minima, which in part due to the high curvature of the mesh. From the oscillation
of the loss in Table 6 there are signs that λ is set too high, but this is somewhat
necessary so as to allow enough movement past local minima. The gradient ascent
therefore seems much more useful in determining the fine alignment, after a cruder
alignment from, for example, the adaptive grid search.

1At least if sampling uniformly at random.

20 HENRY KIRVESLAHTI AND XIAOHAN WANG

7. Code availability
The Ectoplasm code is available on:
https://github.com/hkirvesl/ECTOPLASM.

8. Discussion and Future directions

In this work we have introduced a digitalization of the Euler Characteristic
Transform and developed a functioning algorithm for computing distances with it
in a fully digital manner. This is just a first step to making full use of the rich
mathematical structure in the transform. An interesting an important question is
how one might use this transform for more complicated statistical problems, such
as subshape analysis, in the spirit of [14].
One of the key enablers of the digitalization algorithm is the inclusion exclusion
property of the Euler Characteristic. In theory, it would be quite straightforward to
apply the same procedure for the Persistent Homology Transform. However, com-
puting it in practice is a much harder problem due to the complexity of computa-
tions. But this would be very interesting due to the connections to the sheafification
of shape space as outlined in [1].
The differentiability of the alignment problem provides another connection to
studying the shape space. In this work, we only considered transformations in the
group SO(3) for the specific problem of aligning shapes. But maybe we could find
another more flexible group of transforms that might give us a way to find paths
between shapes via gradients. A simple minded idea would be to differentiate the
penalty with respect to the vertex locations, but care needs to be taken so as not
to allow for self-intersections.
Another problem of applied interest would be extending the algorithm to soft
tissue, as in [8]. This is again theoretically conceivable, but hard in practice.

Acknowledgements
The authors would like to thank Nicolas Berkouk, Sayan Mukherjee, Justin
Curry, Shreya Arya, Kathryn Hess, Robert Ravier, Wojciech Chachólski and Heather
Harrington for helpful discussions and comments.
HK developed a 2-dimensional prototype of the algorithm during his visit to the
Quantum Center at University of Southern Denmark in 2022. The author would
like to thank the Center for their hospitality. The bulk of the work was carried out
at EPFL.

References
[1] Shreya Arya, Justin Curry, and Sayan Mukherjee. A sheaf-theoretic construction of shape
space. Foundations of Computational Mathematics, pages 1–51, 2024.
[2] Mladen Bestvina and Noel Brady. Morse theory and finiteness properties of groups. Inven-
tiones mathematicae, 129:445–470, 1997.
[3] Doug M Boyer, Yaron Lipman, Elizabeth St. Clair, Jesus Puente, Biren A Patel, Thomas
Funkhouser, Jukka Jernvall, and Ingrid Daubechies. Algorithms to automatically quantify the
geometric similarity of anatomical surfaces. Proceedings of the National Academy of Sciences,
108(45):18221–18226, 2011.
[4] Doug M Boyer, Jesus Puente, Justin T Gladman, Chris Glynn, Sayan Mukherjee, Gabriel S
Yapuncich, and Ingrid Daubechies. A new fully automated approach for aligning and com-
paring shapes. The Anatomical Record, 298(1):249–276, 2015.
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 21

[5] Justin Curry, Sayan Mukherjee, and Katharine Turner. How many directions determine a
shape and other sufficiency results for two topological transforms, 05 2018.
[6] Robert Israel (https://mathoverflow.net/users/13650/robert israel). Clairaut’s re-
lation and the equation of great circle in spherical coordinates. MathOverflow.
URL:https://mathoverflow.net/q/247336 (version: 2016-08-12).
[7] Lyman M. Kells, Willis F. Kern, and James R. Bland. Plane And Spherical Trigonometry.
McGraw Hill Book Company, Inc., 1940. Retrieved July 13, 2018.
[8] Henry Kirveslahti and Sayan Mukherjee. Representing fields without correspondences: the
lifted euler characteristic transform. Journal of Applied and Computational Topology, 8(1):1–
34, 2024.
[9] Jacob Leygonie, Steve Oudot, and Ulrike Tillmann. A framework for differential calculus on
persistence barcodes. Foundations of Computational Mathematics, pages 1–63, 2022.
[10] Jesus Puente. Distances and algorithms to compare sets of shapes for automated biological
morphometrics, 2013.
[11] Pierre Schapira. Tomography of constructible functions. In Gérard Cohen, Marc Giusti, and
Teo Mora, editors, Applied Algebra, Algebraic Algorithms and Error-Correcting Codes, pages
427–435, Berlin, Heidelberg, 1995. Springer Berlin Heidelberg.
[12] 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.
[13] Lou van den Dries. Tame Topology and O-minimal Structures. London Mathematical Society
Lecture Note Series. Cambridge University Press, 1998.
[14] Bruce Wang, Timothy Sudijono, Henry Kirveslahti, Tingran Gao, Douglas M Boyer, Sayan
Mukherjee, and Lorin Crawford. A statistical pipeline for identifying physical features that
differentiate classes of 3d shapes. The Annals of Applied Statistics, 15(2):638–661, 2021.

Appendix A. Spherical integration formulae

In this Appendix we detail how we integrate the height functions over spherical
polygons in three dimensions.
We begin by fixing spherical coordinates for points on the unit sphere S 2 , i.e. a
homeomorphism
F : (0, 2π) × (−π, π) → S 2 − I, (ϕ, τ ) 7→ (cos τ cos ϕ, cos τ sin ϕ, sin τ ),
where I is the union of Prime Meridian, North Pole and South Pole. In geometry, ϕ
is the azimuthal angle and τ is the signed angle between the vector from the origin
to that point and the xy-plane. We now provide a parametrization of great circles
under this coordinate system.
Lemma A.1 ([6]). Under the chosen spherical coordinate, a non-polar great circle
on a unit sphere can be characterized by
tan τ = a cos(ϕ − ϕ0 ),
in which we denote a = tan τ0 and (ϕ0 , τ0 ) refers to the point with the minimum
z-coordinate on the great circle.
Proof. Firstly note that tan τ = 0 characterizes the equator, where a = tanτ0 = 0.
If the great circle, denoted by G, is not the equator, there exists only one point
with the minimum z-coordinate on it. Applying a counterclockwise rotation by ϕ0
to G, denoted by ρϕ0 G, we notice that the unit vector u = [sin(−τ0 ), 0, cos(−τ0 )]
is normal to the plane that ρϕ0 G lies on. Therefore
u · [cos τ cos (ϕ − ϕ0 ), cos τ sin (ϕ − ϕ0 ), sin τ ] = 0,
and hence we get tan τ = a cos(ϕ − ϕ0 ). □
22 HENRY KIRVESLAHTI AND XIAOHAN WANG

With this parametrization for great circles, we can compute the following spherical
integration we need for spherical polygons.
Proposition A.2. Given a spherical polygon SPi = {(ϕ1 , τ1 ), (ϕ2 , τ2 ), ..., (ϕn , τn )}
and associated vertex pi = (ϕi , τi ), the ordering of vertices of SPi is compatible with
∂ ∂
the orientation of the polygon, namely [ ∂ϕ , ∂τ ]. If (SPi , pi ) satisfies
1. pi ∈
/ I;
2. SPi does not intersect with I, equator and any meridian.
we denote ϕn+1 = ϕ1 and then have
Z n
X 1 1 1
I(pi ) = opi · v dϕdτ = sin τi · Ik1 − cos τi · Ik2 − cos τi · Ik3 ,
SPi 4 4 2
k=1

where we set the origin to be o, v is the direction vector and

Z ϕk+1
1 − a2 cos2 (ϕ − ϕ0 )
Ik1 = dϕ,
ϕk 1 − a2 cos2 (ϕ − ϕ0 )
Z ϕk+1
2a · cos(ϕ − ϕ0 ) cos(ϕ − ϕi )
Ik2 = dϕ,
ϕk 1 + a2 cos2 (ϕ − ϕ0 )
Z ϕk+1
Ik3 = arctan(a · cos(ϕ − ϕ0 )) cos(ϕ − ϕi ) dϕ.
ϕk

Proof. According to [7], the angle σ between direction vector v = (ϕ, τ ) and opi
satisfies
cos σ = sin τi sin τ + cos τi cos τ cos(ϕ − ϕi ).
Notice that opi · v = cos σ, we have
Z Z
opi · v dϕdτ = sin τi sin τ + cos τi cos τ cos(ϕ − ϕi ) dϕdτ.
SPi SPi

Denote ω = sin τi · cos 2τ − cos τi · sin 2τ · cos(ϕ − ϕi ) − 21 cos τi · cos(ϕ − ϕi ) the

1
4
1
4
1-form on SPi , we observe that dω is exactly the 2-form being integrated on the
right-hand side. By Stokes Theorem, we have
Z Z Z Xn Z
opi · v dϕdτ = dω = ω= ω,
SPi SPi ∂(SPi ) k=1 ek

where ek is the edge of SPi from (ϕk , τk ) to (ϕk+1 , τk+1 ). We then use the tangent
half-angle formula to replace cos 2τ , sin 2τ by tan τ and plug the parametrization
in lemma A.1 into ω.
1 − tan2 τ
Z Z 1 1 2 tan τ 1 
ω= sin τi · 2τ
− cos τi 2τ
cos (ϕ − ϕi ) − cos τi · τ · cos (ϕ − ϕi ) dϕ
ek ek 4 1 + tan 4 1 + tan 2
1 1 1
= sin τi · Ik1 − cos τi · Ik2 − cos τi · Ik3 .
4 4 2
□

As we can see, there are restrictions on the choice of (SPi , pi ) to make the coor-
dinate system and parametrization applicable. In practice, it is always possible
and efficient to apply a rotation to (SPi , pi ) to satisfy the restrictions. Hence the
 DIGITAL TOOLS FOR NON-DIFFEOMORPHIC SHAPES 23

formula above is enough for us to compute the desired integration on any kind of
spherical polygon.
Laboratory for Topology and Neuroscience, EPFL, Lausanne, Switzerland; IMADA
& Center for Quantum Mathematics, University of Southern Denmark, Odense, Den-
mark
Email address: hklahti@imada.sdu.dk

Laboratory for Topology and Neuroscience, EPFL, Lausanne, Switzerland; School

of Physical and Mathematical Sciences, Nanyang Technological University, Singapore
Email address: xiaohan005@e.ntu.edu.sg