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

HAL
open science

# Cellular Fourier analysis for geometrically disordered materials

Antoine Fruleux, Arezki Boudaoud

► To cite this version:

Antoine Fruleux, Arezki Boudaoud. Cellular Fourier analysis for geometrically disordered materials. Physical Review Research, 2021, 3 (2), 10.1103/PhysRevResearch.3.023036. hal-03298520

HAL Id: hal-03298520

https://hal.inrae.fr/hal-03298520v1

Submitted on 20 Jan 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cellular Fourier analysis for geometrically disordered materials

Antoine Fruleux* and Arezki Boudaoud†

RDP, Université de Lyon, ENS de Lyon, UCB Lyon 1, INRAE, CNRS, 69364 Lyon Cedex 07, France and LadHyX, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, 91128 Palaiseau Cedex, France

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

(Received 19 October 2020; accepted 18 March 2021; published 12 April 2021)

Many media are divided into elementary units with irregular shape and size, as exemplified by domains in magnetic materials, bubbles in foams, or cells in biological tissues. Such media are essentially characterized by geometrical disorder of their elementary units, which we term cells. Cells set a reference scale at which are often assessed parameters and fields reflecting material properties and state. Here, we consider the spectral analysis of spatially varying fields. Such analysis is difficult in geometrically disordered media, because space discretization based on standard coordinate systems is not commensurate with the natural discretization into geometrically disordered cells. Indeed, we found that two classical spectral methods, the fast Fourier transform and the graph Fourier transform, fail to reproduce all expected properties of spectra of plane waves and of white noise. We therefore built a method, which we call cellular Fourier transform (CFT), to analyze cell-scale fields, which comprise both discrete fields defined only at cell level and continuous fields smoothed out from their subcell variations. Our approach is based on the construction of a discrete operator suited to the disordered geometry and on the computation of its eigenvectors, which, respectively, play the same role as the Laplace operator and sine waves in Euclidean coordinate systems. We show that CFT has the expected behavior for sinusoidal fields and for random fields with long-range correlations. Our approach for spectral analysis is suited to any geometrically disordered material, such as a biological tissue with complex geometry, opening the path to systematic multiscale analyses of material behavior.

DOI: 10.1103/PhysRevResearch.3.023036

# I. INTRODUCTION

The past decades have seen a growing interest in geometrically disordered media [1] such as liquid and solid foams [2–4], granular materials [5], or biological tissues [6,7]. This brought many questions and concepts related to the dynamics of these media such as coarsening [8–13], fluctuations [14], jamming transition [15–19], grain growth [20], or applicability to living tissues [21]. Many experimental approaches were developed to observe and quantify cell tilings in these media. For instance, magnetic resonance imaging [22] or x-ray tomography [11,23] enable imaging of foam evolution in three dimensions. Imaging of biological tissue is performed with serial block-face scanning electron microscopy [24] or with confocal microscopy of living samples [21]. Using efficient algorithms such as the watershed transform [25], it has been possible to segment these two-dimensional (2D) and three-dimensional (3D) images, i.e., to extract the geometry and the arrangement of the cells, as performed in foams [26], granular material [27], or in biological tissues [28]. Here, we consider

quantitative analyses of properties or fields defined on such segmented images.

An example is given in Fig. 1, which shows a coarsening 2D liquid foam. It is disordered, constituted of cells of broadly distributed sizes, and irregularly arranged. Accordingly, analyses of fields defined in geometrically disordered media require disentangling potential randomness associated with the field from randomness due to geometry. In addition, the cell often provides a reference scale below which it is difficult, irrelevant, or impossible to define fields. In a foam, areal growth is defined at a discrete level (generally at cell level, Fig. 1) because growth requires landmarks (here, vertices) to be computed. These special features make it difficult to assess spatial patterns and test theoretical predictions based on continuous models, such as our prediction of long-range spatial correlations for growth fluctuations in biological tissues [29]. Here, we develop an approach to overcome these difficulties, based on harmonic representation of signals defined on cellular media, which enables one to properly analyze the spectra of these signals. We term our approach cellular Fourier transform (CFT).

Spectral analysis decomposes signals into linear combinations of harmonics [30]. Ad hoc harmonics depend on how signals are represented [31]: for a continuous signal in Euclidean space, it is common to use plane waves and to consider the Fourier transform, while discrete signals defined on regular grids are often decomposed into the eigenvectors of circulant matrices, yielding the fast Fourier transform (FFT). Initial frameworks for spectral analysis of signals on

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>* antoine.fruleux@polytechnique.edu</footnote>

<footnote>† arezki.boudaoud@polytechnique.edu</footnote>
Areal growth rate in 2D foam

<figure><img src="image_4.png" /><figcaption>FIG. 1. Cell areal growth rate in a 2D coarsening foam with visible polydispersity in size. Black lines represent the liquid films between cells, which are colored according to their relative areal growth rate computed over a period of 106 s (color scale on right). Data courtesy of Jérôme Duplat (see [13]).</figcaption></figure>

irregular grids were based on FFT [32]. More recent approaches constructed ad hoc harmonics on graphs [33]. Harmonics depend on geometry and, in compact metric spaces, they can be defined as eigenfunctions of the Laplace operator [30]. This idea has been extensively used in discrete analysis [34], especially to analyze signals on graphs [33,35,36]. A graph may be endowed with an irregular geometry by ascribing a distance to each edge and define Laplace operators that incorporate distances. However, graphs cannot be used to describe geometrically disordered materials, because graphs do not account for the full geometry of unit cells. Discrete Laplace operators have also been defined for triangular meshing of surfaces [37], but their use for signals defined on geometrically disordered materials seems problematic for several reasons. These operators converge only weakly to the smooth Laplace-Beltrami operator in the limit of small mesh size [38] and discrete Laplacians on triangular meshes cannot satisfy all desired natural properties [39]. Other reasons are related to the nature of the mathematical object that we consider: Cellular tessellations of space (2D or 3D), in which cells may have complex shapes and varying topologies (number of neighbors). Therefore the arrangement of cells cannot be encapsulated in binary relations as in triangular meshes or in graphs.

We therefore built a framework to analyze signals defined on (possibly disordered) tessellations of space. We first motivate our study by applying the fast Fourier transform and the graph Fourier transform to examples of cellularized signals. We present the geometry of the medium and how signals are represented. We define a coarse Laplace operator, applicable to signals with variations at sub- and supracell scale; we show that sine waves are its eigenfunction in the Euclidean space. We project this operator on the cellularized geometry and discretize it. Finally, we test our analysis with numerically generated data, illustrate it with experimental data from a coarsening foam, and discuss the potential applications of this framework. For simplicity, we present results for polygonal tilings of the Euclidean plane, but this method is broadly applicable to domains of any geometry and dimension. We also provide MATLAB (Mathworks) scripts implementing CFT

that can be readily adapted to analyze any type of cellularized signal.

II. MOTIVATION

We first discuss the applicability of the FFT or the GFT to
cellularized signals.

A. Reference tessellation and reference fields

To test the different approaches, we generated a cellular tessellation and reference fields that associate to each cell a scalar value. We started from 1000 × 1000 images of a range of continuous fields shown in Figs. 2(a)–2(d): plane waves of fixed amplitude and four values of wave number, step function, long-range correlated noise, and white noise. Stationary waves [Fig. 2(a)] have wave numbers $2\sqrt{2\pi n}/a$ with $n = 1, 3, 5,$ and 10. The step function is defined as $f(x, y) = -1$ if $y < 0$ and $f(x, y) = 1$ if $y \ge 0$ [Fig. 2(b)]. We built the long-range correlated random field [Fig. 2(c)] using the Fourier filtering method detailed in [40]; the two-point correlation function of pixel intensity $\langle I_i I_j \rangle$ decays with the distance $d_{ij}$ between pixels like $1/d_{ij}$. The white noise has amplitude 1 [Fig. 2(d)].

We built a reference domain and its partition into cells by
generating the Voronoi tessellation of a random distribution
of 1000 seeds initially placed in a square of side 1000. We
only kept the 894 cells entirely included in the square, yielding
the domain Ω as shown in Figs. 2(e)–2(h). Ω is not a perfect
square; it is a polygon with a large number of edges. The
number of cell neighbors is broadly distributed around six. For
a given continuous field, we computed the average of the field
over each cell; the mapping between a cell and the average
field over that cell defines the cellularized field. The four types
of cellularized fields are shown in Figs. 2(e)–2(h) and appear
as piecewise constant versions of the continuous fields. In the
remainder of this section, we focus on stationary waves and on
white noise [Figs. 2(a), 2(d) 2(e), and 2(h)]; we use the step
function and the long-range correlated noise in Sec. IV.

B. Applicability of fast Fourier transform

We first applied the FFT to reference fields. We discretized
these fields on 1000 × 1000 images, setting to zero the in-
tensity of pixels not belonging to cells (due to the Voronoi
tessellation not filling the square). We show the spectra of the
reference fields as a function of the x and y components of the
wave vector [Fig. 3(a) for stationary waves and Fig. 3(b) for
white noise] and as a function of the wave number in Fig. 3(c).
We made wave vectors dimensionless using the typical cell
scale $l_c$, defined as the square root of the mean cell area.
To ease readability of figures, we only plotted the first few
thousands of modes among the million of possible modes (the
starting data has dimensions 1000 × 1000).

The spectra of the stationary waves show four peaks
(±q<sub>x</sub><sup>(m)</sup>, ±q<sub>y</sub><sup>(m)</sup>), lines joining them, and background noise
[Fig. 3(a)]. The four peaks correspond to the wave numbers
of the initial stationary wave while the lines result from a
windowing effect due to the Voronoi tessellation not entirely
filling the square. The spectrum of the white noise is maximal
for the wave vector (0,0) and decays slowly away from (0,0)
<figure><img src="image_5.png" /><figcaption>FIG. 2. Reference fields: four types of continuous fields and their cellularized versions defined on a reference domain. (a), (b), (c), and (d) respectively show stationary waves, a step function, a random field with long-range correlation, and a white noise, while (e), (f), (g), and (h) show their respective cellularized versions, in which each cell is ascribed the average of the continuous field over that cell.</figcaption></figure>

[Fig. 3(b)]. In all cases, the amplitude of the FFT increasingly departs from expectations as the wave number increases: the height of the peaks decrease for waves of increasing wave numbers [in Fig. 3(c)], instead of having a fixed amplitude, and the spectrum of the white noise decreases with wave number [in Fig. 3(d)] instead of having a constant amplitude, a phenomenon that is very clear when averaging over realizations of the noise [inset of Fig. 3(d)].

The higher the wave number, the less the signal can be approximated by its cellularized projection. Cells act like a low-pass filter that reduces the amplitude of higher FFT modes with respect to the spectra of the initial signals. As we will see later, the method that we developed circumvents this issue.

### C. Applicability of graph Fourier transform

The GFT is applicable to signals that are defined on graphs; the harmonics are then the eigenmodes of a discrete Laplace operator. Here a graph can be defined from cells and neighborhood relations between cells: each vertex of the graph corresponds to one cell, and this vertex is connected to all vertices corresponding to neighboring cells. We may account for the metrics of the cellular tessellation by associating to each edge (link between graph vertices *i* and *j*) a weight $w_{ij}$ that depends on the distance between cell centers $d_{ij}$, $w_{ij} = f(d_{ij})$; the weight $w_{ij}$ vanishes if the two cells are not neighbors. For a field $\phi$ that has values $\phi_i$ for *i* spanning vertices, the value of the Laplacian at vertex *i* is $[\mathcal{L}(\phi)]_i = \sum_j w_{ij}(\phi_j - \phi)$. For a fair comparison with the method that we developed, we used the same kernel $f(r) = \exp(-r/\sigma)$ with $\sigma = 7l_c$. We computed the eigenvectors $e_k$ of the Laplacian $\mathcal{L}$, which are associated to the eigenvalues $\hat{L}_k$. The graph Fourier transform of a field is then given by the components of the field with respect to the basis $\{e_k\}$.

Because there is no standard definition of the wave number for GFT harmonics, we plotted the spectra as a function of the eigenvalue $\hat{L}_k$ of the graph Laplacian $\mathcal{L}$ [Figs. 3(e)–3(f)]. The GFT behaves poorly for stationary waves, having a marked peak only for the wave with highest wavelength. The GFT spectrum of the white noise has amplitudes that are one order of magnitude lower than expected; in addition, when the

spectrum is averaged over realizations of the noise, a slowly decreasing trend appears in the spectrum, in contradiction with expected constant amplitude.

The GFT is not really applicable to cellularized fields because it behaves poorly with respect to spectra of the continuous fields and it is difficult to define wave numbers. An explanation is that geometrical information on cell shape is partially lost when considering a weight that depends on distance between cell centers. As we will see later, the method that we developed fully accounts for cell geometry.

## III. FORMULATION

We first present the mathematical basis of the CFT. The reader may skip to Sec. IV A to find the implementation and validation of the CFT. The main results are given by Eqs. (13)–(16) that, together with Eq. (4), define the harmonics and their wave numbers from an appropriate discrete Laplace operator.

### A. Signal representation on a cellularized space

In this section, we explain how signals are represented in the cellularized space and we specify desired properties of the harmonics.

We consider a bounded domain $\Omega$ of the $n$-dimensional space $\mathbb{R}^n$, divided into $N$ subdomains $\{\omega_i, i = 0, 1, \dots, N-1\}$, which we call cells. We consider an integration measure $d\mu$ defined on $\mathbb{R}^n$. For a Euclidean metric and coordinates $(y_1, y_2, \dots, y_n)$, we would have $d\mu(y) = dy_1dy_2\dots dy_n$. Let $S$ and $S_i$ be the sizes (e.g., areas) of domain $\Omega$ and cells $\omega_i$, respectively, $S = \int_{\Omega} 1 d\mu(x)$ and $S_i = \int_{\omega_i} 1 d\mu(x)$, where the integrals over $\Omega$ and the cell domains $\omega_i$ are computed according to the measure $d\mu(y)$. We aim at analyzing a signal $f$ defined over $\Omega$ (and assumed to have sufficient regularity for all formulations to be mathematically well posed). The smoothed version of the signal $f$ is given by a cellwise constant function $f_C$. $f_C$ associates each cell $\omega_i$ with the average of the continuous signal $f$ over that cell, $1/S_i \int_{\omega_i} f(y)d\mu(y)$. We call $f_C$ the representation of $f$. Figure 2 shows examples of signals $f$ and their representations $f_C$ on a random Voronoi tessellation. The representation $f_C$ belongs to a vectorial space
<figure><img src="image_6.png" /><figcaption>FIG. 3. Fast Fourier transform and graph Fourier transform applied to cellular representations of plane waves and white noise [Figs. 2(e) and 2(h)]. FFT modulus is plotted as a function of the dimensionless components of wave vector ($q_x, q_y$) in (a) and (b) and as a function of wave number $q = \sqrt{q_x^2 + q_y^2}$ in (c) and (d), using typical cell scale $l_c$ as unit length. In (e) and (f), GFT modulus is plotted as a function associated eigenvalue $\hat{L}_k$ of Laplacian. Insets in (d) and (f) show the root mean square spectra obtained by averaging 1000 realizations of the white noise. Dashed vertical lines in (c) are at the wave numbers of the initial fields, while dark horizontal lines in (d) and (f) correspond to the spectral density of the initial noise (theoretical curve).</figcaption></figure>

$\mathcal{E}$ of dimension $N$ (the number of cells $\omega_i$), which we call the representation space. A basis of $\mathcal{E}$ is the set of functions $\{\psi_i\}$ that vanish outside $\omega_i$ and are defined by

$$ \psi_i: \begin{array}{l} \Omega \to \mathbb{R}, \\ x \to \psi_i(x) = \begin{cases} 1/\sqrt{S_i} & \text{if } x \in \omega_i \\ 0 & \text{else.} \end{cases} \end{array} \quad (1) $$

Given the standard scalar product $\langle f \cdot g \rangle = \int_{\Omega} f(x)g(x)d\mu(x)$ of two functions $f$ and $g$, the basis of functions $\{\psi_i\}$ is orthonormal, $\langle \psi_i \cdot \psi_j \rangle = \delta_{ij}$, where $\delta_{ij}$ is the Kronecker symbol ($\delta_{ij} = 1$ if $i = j$, else $\delta_{ij} = 0$). The representation $f_C$ of $f$ is also its orthogonal projection on the representation space $\mathcal{E}$,

$$ f_C(x) = \sum_{i=0}^{N-1} f_i \psi_i(x) \quad \text{with} \quad f_i = \langle \psi_i \cdot f \rangle. \quad (2) $$

We note that $f_i = 1/\sqrt{S_i} \int_{\omega_i} d\mu(y)f(y)$ is not the signal averaged over the domain $\omega_i$. The prefactor involving size $S_i$ of cell $i$ is important for the present analysis to be mathematically sound and the definition of $f_i$ is the appropriate discretization of the cellularized signal. Inappropriate discretization partly explains the poor results of GFT. In the following, we aim to define another orthogonal basis for the representation space $\mathcal{E}$ so that its elements enable spectral analysis in cellular media. We will call the elements $e_k$ of this basis harmonics of the representation space. They can be written as $e_k(x) = \sum_i U_{ki} \psi_i(x)$, where $U_{ki} = \langle e_k \cdot \psi_i \rangle$ are the elements of the rotation matrix $U$ between the two bases. Finding the harmonics of the representation space is equivalent to determining the unitary matrix $U$.

## B. Coarse Laplace operator

In infinite Euclidean space, Fourier harmonics are plane waves. These plane waves are notably eigenfunctions of the classical Laplace operator. More generally, these plane waves are eigenfunctions of all integral operators that are invariant by translation, a property that we will use to define the harmonics $e_k$. In this section we consider a Laplace-like integral operator and investigate its properties, first in infinite Euclidean space and then in bounded domain; we then explain how the problem can be discretized to define the harmonics $e_k$.

### 1. In unbounded space

Harmonics are often defined as eigenfunctions of the Laplace operator. Because the signals that we consider are smoothed out of their subcellular variations, we build a coarse version of the Laplace operator, formulated as an integral operator $\mathcal{L}$ which, to each function $f$ defined on $\mathbb{R}^n$, associates

$$ \mathcal{L}[f](x) = \int_{\mathbb{R}^n} w(|x-z|)[f(x) - f(z)]d\mu(z). \quad (3) $$

The kernel $w$ is an integrable real function and $|x - z|$ is the Euclidean distance between points $x$ and $z$ of $\mathbb{R}^n$. We assume that $w(r)$ has a maximum at $r = 0$ and vanishes when $r \to \infty$, with a characteristic decay length $\sigma$. Like the discrete Laplacian on a grid, the operator $\mathcal{W}$ averages the difference between the local field $f(x)$ and its value on the neighborhood of $x$. In the limit where the length scale $\sigma$ vanishes, $\mathcal{L}[f] \simeq C\nabla^2 f$, where $\nabla^2$ is the classical Laplace operator and the constant $C = -\int_{\mathbb{R}^n} w(|z|)z^2/2 d\mu(z)$. The operator $\mathcal{L}$ can be seen as a coarse version of the Laplacian $\nabla^2$ (see [41]).

In all generality, plane waves of wave number $q$, $u_q(x) = \exp(Iq \cdot x)$ (with $I^2 = -1$) are eigenfunctions of $\mathcal{L}$:

$$ \begin{align} \mathcal{L}[u_q](x) &= \hat{L}(|q|)u_q(x), \quad \text{with} \nonumber \\ \hat{L}(|q|) &= \hat{W}(0) - \hat{W}(|q|), \nonumber \\ \hat{W}(|q|) &= \int_0^{+\infty} w(r)A_d(|q|r)r^{d-1}dr, \quad \text{and} \nonumber \\ A_n(r) &= \frac{2\pi^{(n-1)/2}}{\Gamma((n-1)/2)} \int_0^{\pi} d\theta(\sin\theta)^{n-2} \exp(I|q|r\cos\theta). \tag{4} \end{align} $$
We will later consider the case when $\hat{L}$ is positive and monotonously increasing on $[0, +\infty[$ so that the eigenfunctions of $\mathcal{L}$ associated to a given eigenvalue are linear combinations of plane waves with all wave numbers having the norm $|q|$.

## 2. In bounded domain

We now consider a compact bounded domain $\Omega$ of $\mathbb{R}^n$ and functions $f$ defined on $\Omega$, and we replace the operator of the previous section by

$$ \mathcal{L}[f](x) = \int_{\Omega} w(|x - y|)[f(x) - f(y)]d\mu(y), \quad (5) $$

where $\mathcal{L}$ is a compact and self-adjoint operator; according to the spectral theorem, the eigenfunctions of $\mathcal{L}$ form an orthogonal basis of the associated Hilbert space [30]. If $w$ is well localized around 0 (equivalently if its decay length $\sigma$ is small with respect to size of the domain $\Omega$), then the integral in (5) can be approximated by an integral over the whole Euclidean space $\mathbb{R}^n$ provided that the position $x$ is not too close to the boundary of $\Omega$. Therefore, Eq. (4) implies that the eigenvalues of $\mathcal{L}$ are approximately $\hat{L}(|q|)$. Associated eigenfunctions are locally approximated by linear combinations of plane waves with wave numbers of norm $|q|$, except close to the boundary; eigenfunctions may show boundary layers of width $\sigma$ (see below). Note that, in order to reduce this boundary effect for eigenfunctions associated to the lowest eigenvalues, we defined $\mathcal{L}[f](x)$ in Eq. (5) using the difference $f(x) - f(y)$, which enforces constant functions to be eigenfunctions of $\mathcal{L}$ that are associated to the eigenvalue 0.

## C. Harmonics of the representation space

In this section we discretize the problem and consider the representation of the coarse Laplace operator, $\mathcal{L}$, defined above. We introduce in Sec. III C 1 the harmonics of the representation space and investigate in Sec. III C 2 how they are related to the eigenfunctions of $\mathcal{L}$ and how their wave numbers can be computed. In Sec. III C 3, we introduce a correction to reduce boundary effects.

### 1. Discrete Laplace operator

In order to build the representation of the coarse Laplace operator, we first consider the representation $w_C(x, y)$ of the kernel $w(x, y)$ assumed to be a function of the distance $|x - y|$ so that $w(x, y) = w(|x - y|)$. $w(x, y)$ is a function of two variables and its representation is piecewise constant on $\Omega^2$. Accordingly, we generalize the representation of a function of a single variable defined as a projection on $\mathcal{E}$ and define the projection of $w(x, y)$ on $\mathcal{E}^2$ as follows:

$$ w_C(x, y) = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} W_{ij} \psi_i(x) \psi_j(y), \quad (6) $$

where the elements of the weight matrix $W$ are given by

$$ W_{ij} = \frac{1}{\sqrt{S_i S_j}} \int_{\omega_i} \int_{\omega_j} w(|x - y|) d\mu(x) d\mu(y). \quad (7) $$

The integral operator $\mathcal{L}_C$ associated to the kernel $w_C$ is $\mathcal{L}_C[f](x) = \int_{\Omega} w_C(x, y)[f(x) - f(y)]d\mu(y)$. Using

Eq. (6) for $w_C$ and the relations $\int_{\Omega} \psi_j(y)d\mu(y) = \sqrt{S_j}$ and $f_j = \langle \psi_j \cdot f \rangle$, we rewrite the integral operator as $\mathcal{L}_C[f](x) = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \psi_i(x)W_{ij}[f(x)\sqrt{S_j} - f_j]$. We then apply this operator to a function $f_C(x) = \sum_{i=0}^{N-1} f_i \psi_i(x)$ of the representation space and use the relation $\psi_i(x)\psi_j(x) = (1/\sqrt{S_i})\delta_{ij}\psi_i(x)$, so as to get $\mathcal{L}_C[f_C](x) = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} L_{ij} f_j \psi_i(x)$, where

$$ L = D - W, \quad D_{ij} = \delta_{ij} \sum_{k=0}^{N-1} W_{ik} \sqrt{\frac{S_k}{S_i}}. \quad (8) $$

We call the matrix $L$ the discrete Laplace operator. It shows similarities with classical discrete Laplace operators [39], such as being symmetric, having positive weights $W_{ij}$, and being positive semidefinite. We nevertheless emphasize that these classical discrete Laplacians do not operate on the same mathematical objects: they apply to functions defined on meshes or on graphs, whereas our operator $L$ applies to functions that are piecewise constant over each cell.

### 2. Eigenfunctions and eigenvalues

The discrete Laplace operator $L$ can be diagonalized and we propose to define the unitary matrix $U$ introduced in Sec. III A from the eigenvectors of $L$, yielding

$$ L_{ij} = \sum_{k=0}^{N-1} \hat{L}_k U_{ki} U_{kj}, \quad (9) $$

where $\{\hat{L}_k, k = 0, 1, \dots, N-1\}$ are the eigenvalues of $L$. Note that the eigenvector associated to eigenvalue 0 is $(\sqrt{S_0}, \sqrt{S_1}, \dots, \sqrt{S_{N-1}})$ instead of being $(1, 1, \dots, 1)$ as for classical discrete Laplace operators. This difference illustrates how cell size heterogeneity is taken into account by the cellular Fourier transform. The harmonics $\{e_k, k = 0, 1, \dots, N-1\}$ associated to $U$ are the eigenfunctions of $\mathcal{L}_C$:

$$ \mathcal{L}_C[e_k](x) = \hat{L}_k e_k(x). \quad (10) $$

In the following, we show that harmonics of the representation space are good representations of harmonics of the whole space of functions defined on $\Omega$ and we determine the wave numbers associated with the eigenfunctions $e_k$, by using the coarse Laplace operator $\mathcal{L}$. To do so, we first consider an eigenfunction $f$ of $\mathcal{L}: \mathcal{L}[f](x) = \lambda f(x)$. Based on the definitions above, its representation $f_C$ verifies

$$ \mathcal{L}_C[f_C](x) = \lambda f_C(x) + \sum_{i=0}^{N-1} \frac{\psi_i(x)}{\sqrt{S_i}} \int_{\Omega_i} \int_{\Omega} d\mu(y)w(|z - y|) \times [f(y) - f_C(y)]d\mu(z). \quad (11) $$

In the integral in (11), $f(y) - f_C(y)$ varies quickly, i.e., at the typical cell scale $l_c = (\mu/N)^{1/n}$. If the decay length $\sigma$ of the kernel $w$ is greater than or comparable to $l_c$, then the integral in (11), yielding $\mathcal{L}_C[f_C](x) = \lambda f_C(x)$, and $f_C$ is, within a good approximation, an eigenfunction of $\mathcal{L}_C$, associated to the same eigenvalue $\lambda$. Therefore the linear space generated by the eigenfunctions of $\mathcal{L}$ associated with eigenvalue $\hat{L}(|q_k|)$ [see Eq. (4)] is projected to the linear space generated by the
harmonics $e_k$ such that

$$ \hat{L}_k = \hat{L}(q_k). \qquad (12) $$

Finally, based on Sec. III B, each eigenfunction of $\mathcal{L}$ associated to eigenvalue $\hat{L}(|q_k|)$ is locally well approximated by a linear combination of plane waves with wave vectors having the same norm $|q_k|$. Therefore the $e_k$ are appropriate harmonics of the representation space because they are the projections of linear combinations of plane waves of the same wave number $q_k$, defined by Eq. (12).

### 3. Correcting boundary effects

The results obtained in Secs. III C 1 and 2 are only valid in the bulk of domain $\Omega$. Indeed, when we explicitly compute the set $\{e_k, k = 0, 1, \dots, N-1\}$ following the above, we found that some of the harmonics $e_k$ show visible boundary layers near the edges of $\Omega$. This can be explained qualitatively as follows: At edges, $\mathcal{L}[f](x)$ is an integral over a domain about twice smaller than in the bulk; consequently, variations of eigenfunctions at edges are about twice as big as in bulk. Two qualitatively similar methods could be used to correct this artifact. One is to rescale the weight function $w(|x-y|)$ by its integral over $\Omega$, i.e., use the kernel $w(|x-y|)/\int_{\Omega} w(|x-z|)d\mu(z)$. Another is to directly rescale each row $i$ of the discrete Laplace operator by $D_{ii} = \frac{1}{S_i} \int_{\omega_i} \int_{\Omega} w(|x-y|)d\mu(y)d\mu(x)$. Here we use the latter method and we consider the rescaled Laplace $\bar{L}$ operator,

$$ \bar{L}_{ij} = \delta_{ij} - \bar{W}_{ij}, \quad \text{with} \\ \bar{W}_{ij} = \sqrt{\frac{S_i}{S_j} \frac{\int_{\omega_i} \int_{\omega_j} w(|x-y|)d\mu(x)d\mu(y)}{\int_{\omega_i} \int_{\Omega} w(|z-t|)d\mu(z)d\mu(t)}}. \qquad (13) $$

This rescaling only marginally changes the matrix $L$, though it breaks its symmetry. For this reason $U$ is no longer defined as the eigenvectors of $\bar{L}$ but as its right-singular vectors. This does not significantly affect $U$, except for columns corresponding to cells at the edges of $\Omega$. The singular values $\hat{L}_k$ of the rescaled Laplace operator are therefore still related to wave numbers via Eq. (12). If $V$ are the left-singular vectors, we may write the components of $\bar{L}$ as

$$ \bar{L}_{ij} = \sum_{k=0}^{N-1} \hat{L}_k V_{ki} U_{kj}. \qquad (14) $$

The wave number $q_k$ can be defined from the eigenvalue $\hat{L}_k$ using Eq. (4),

$$ \hat{L}_k = \hat{L}(q_k)/\hat{L}(0). \qquad (15) $$

From this analysis, we can deduce the spatial spectrum of a signal $f$. The $k$th spectral coefficient is given by

$$ \hat{\Phi}_k = \langle e_k \cdot f \rangle = \sum_{i=0}^{N-1} U_{ki} \frac{1}{\sqrt{S_i}} \int_{\omega_i} f(x) d\mu(x), \qquad (16) $$

and is associated to the wave number $q_k$ defined in Eq. (15). We call this spectrum the cellular Fourier transform: It is the appropriate equivalent of the classical Fourier spectrum in $\mathbb{R}^n$ for signals defined at cell level. The spectrum is fully defined by Eqs. (13)–(16), together with Eq. (4).

<figure><img src="image_7.png" /><figcaption>FIG. 4. The relation between wave numbers and eigenvalues of $L$. The black line shows the function $Q$ and the blue crosses show $Q(\hat{L}_k)$ that correspond to the discrete wave numbers. Kernel: exponential kernel, $\sigma = 7l_c$.</figcaption></figure>

## IV. IMPLEMENTATION AND RESULTS

### A. Implementation

We implemented the CFT on the domain and for the fields shown in Fig. 2. The CFT relies on defining harmonics as the eigenvectors of the discrete Laplace operator given by Eq. (13). This equation involves cell sizes $S_i$ (e.g., cell areas in 2D) and a kernel $w(r)$ that vanishes for large $r$. We tested several kernels $w$ and several values of their characteristic decay lengths $\sigma$ (see below). We numerically computed the integrals in Eq. (13) using Gauss quadrature. For tractability of numerical calculations, we used the approximation $w(r > 6\sigma) = 0$. We obtained the harmonics and associated wave numbers [Eqs. (14) and (15)] using the singular value decomposition algorithm of MATLAB (Mathworks) (see Appendix). We present in Sec. IV B the CFTs of a few results of our approach before testing in Sec. IV D how the parameters of the kernel influence the analysis. As an illustration, we also performed an analysis of the foam coarsening data shown in Fig. 1, which is presented in Sec. IV E.

### B. CFT applied to artificial fields

In this section, we use an exponential kernel, $w(x) = \sigma^{-2}\exp(-r/\sigma)$. Following (12), the corresponding wave numbers $q_k$ take values $q_k = 1/\sigma Q(\hat{L}_k)$, where $Q(l) = \sqrt{(1-l)^{-2/3} - 1}$ is the inverse function of $\hat{L}$ introduced in (4). Note that the wave number $q_k$ associated to the harmonics $e_k$ is not given by the square root of $\hat{L}_k$, as would be for the Fourier transform in infinite space (though $q_k \sim \sqrt{\hat{L}_k}$ for small $\hat{L}_k$). These wave numbers are shown in Fig. 4.

Figure 5 shows the first harmonics $e_k$ on $\Omega$. These harmonics resemble linear combinations of the eigenfunctions of the Laplace operator with Neumann boundary conditions on a square of side $a$ [42], $u_{mn}(x, y) = \cos(\pi mx/a) \cos(\pi ny/a)$, where $(x, y)$ are the Cartesian coordinates and $(n, m)$ are integers. Indeed, the shape of $\Omega$ is close to a square and in the limit $\sigma \to 0$, the coarse Laplace operator $\mathcal{L}$ defined in
<figure><img src="image_8.png" /><figcaption>FIG. 5. The six first harmonics on a domain Ω that was generated using a randomly seeded Voronoi tessellation. Blue and yellow correspond, respectively, to negative and positive values of the harmonic. Kernel: exponential, $\sigma = 7l_c$.</figcaption></figure>

Eq. (5) converges toward the classical Laplace operator with Neumann conditions on the boundaries of its domain Ω [41]. Accordingly, the first harmonic of the discrete Laplace operator, $e_0$, is constant like $u_{00}$. The following harmonics $e_1$ and $e_2$ correspond, respectively, to $u_{10}$ and $u_{01}$, but their orientation deviates slightly toward diagonals. More generally, a given harmonic $e_k$ corresponds to a linear combination of eigenfunctions $u_{mn}$ with indices corresponding to $\pi/a\sqrt{m^2+n^2} \sim q_k$. This linear combination is such that higher harmonics $e_k$ tend to have the same spatial periodicity in all directions. Accordingly, the wave numbers of the CFT modes are well defined, but the directions of the wave vectors of the modes are not. Indeed, the CFT was designed for generic geometries, which do not have well-defined wave vectors; for instance, classical harmonics in a disk are Bessel functions of the radial coordinate multiplied by sines of the orthoradial coordinate.

In Fig. 6 we show the spectra of the fields plotted in Fig. 2. As expected, the spectra of the stationary waves representations are maximal close to the wave numbers of the initial continuous waves [Fig. 6(a)]. Nevertheless, the width of the peaks increases with the wave number of the initial continuous stationary waves, owing to a less accurate approximation of the stationary wave by its representation in the face of finite cell size. Also as expected, the spectrum of the step function representation is peaked at 0 [Fig. 6(b)]. The spectra of field with long-range correlations [Fig. 6(c)] and of the white noise [Fig. 6(d)] are random and the amplitudes are distributed around zero (here we plot the absolute value); the amplitudes slowly decay with $q_k$ in C but stay constant in D. To further study random fields, we generated 1000 realizations of the noise. We estimated the average spectrum, which should give an estimate of the Fourier transform of the correlation function. The spectra of representations behave as expected: constant for white noise and power-law decay for field with long-range correlations [Fig. 6(e)].

<figure><img src="image_9.png" /><figcaption>FIG. 6. CFT of the cellularized fields shown in Figs. 2(a)–2(d) Spectra of the cellularized fields shown in Figs. 2(e)–2(h), respectively. (e) Root-mean-square spectra of white noise (red) and of long-range correlated noise (blue); the average is over 1000 realizations and is shown in linear (main plot) and logarithmic (inset) scales. The magenta, yellow, red, and blue data in (a) are associated to the wave numbers $2\sqrt{2\pi n}/a$ of the initial continuous signals (represented by dashed lines) with $n = 1, 3, 5, \text{and } 10$, respectively, and $a$ is the dimension of the square. The dark lines in (c)–(e) represent the theoretical curves given by the Fourier transforms of the continuous signals $f$. Kernel: exponential, $\sigma = 7l_c$.</figcaption></figure>

### C. CFT compared to FFT and GFT

In the context of cellularized signals, the CFT performs better than the fast Fourier transform. They both retrieve peaks in the spectra, but the position of the peaks is less precise with FFT and the peaks are broader for FTT. This can be quantified using the spectral density $\hat{\Phi}_k$; the mean, $q_M$, and root mean square, $q_{rms}$, of wave number are computed as $(\sum_{k=0}^{N-1} q_k^n |\hat{\Phi}_k|^2) / (\sum_{k=0}^{N-1} |\hat{\Phi}_k|^2)$, with $n = 1$ and $n = 2$, respectively. We computed the relative shift in wave number $\delta = (q_M - q_i)/q_i$ between the resulting spectrum and the initial stationary wave. For the CFT, we found values $\delta = -0.0583, -0.127, -0.123$, and $-0.0993$ (for the four respective stationary waves), which are orders of magnitude smaller than for the CFT, $\delta = 3.91, 1.89, 1.89$, and $1.76$. Similarly, we computed the relative standard deviation $\nu = (q_{rms}^2 - q_M^2)^{1/2}/q_i$ of the resulting spectrum. For the CFT, we found values $\nu = 0.164, 0.130, 0.130$, and $0.183$, which are also orders of magnitude lower than for the FFT, $\nu = 26.6, 10.4, 8.09$, and $5.75$. The CFT retrieves a flat spectrum for white noise, whereas the FFT
yields a decaying spectrum. This is due to the FFT adding artificially a large number of degrees of freedom (~$10^6$) to those of the cellularized signal (number of cells ~$10^3$); hence the spectrum of the noise spreads out to larger wave numbers. Reducing the number of FFT modes would not be a solution because FFT is defined on a square grid which cannot match the topological relations between cells in a disordered medium.

In this context, the CFT also performs better than the graph Fourier transform. The GFT yields broad peaks for the spectrum of stationary waves and a decaying spectrum for the white noise. Although this might be improved by a better choice of the weight of links between cells, the GFT would still discard cell shape, and the difficulty of defining wave numbers in the GFT would remain a major issue.

We note, however, that all transforms would yield very similar results for a tessellation with identical square-shaped cells. As such, the CFT is most appropriate for geometrically disordered media.

### D. Sensitivity of the CFT to the kernel

We first tested the influence of the kernel decay length, $\sigma$, of the exponential kernel. We show in Fig. 7 spectra obtained with different values of $\sigma$. We see that when $\sigma$ is small, spectra are shifted toward higher wave numbers. Optimizing the kernel requires taking into account domain size and numerical precision at which the singular value decomposition is performed. The typical cell scale $l_c$ should be a lower bound for $\sigma$ while domain size and numerical precision prevent $\sigma$ from being too large. Optimal values for $\sigma$ must depend on the whole distribution of cell shapes and size as well as on the number of cells. We estimated this optimal value by maximizing the agreement between the wave numbers corresponding to maxima of spectra in Fig. 7(a) and the wave number of the corresponding stationary waves. We found $\sigma \simeq 7l_c$, which is intermediate between cell size and domain size.

We calculated the harmonics of the same domain $\Omega$ using a Gaussian kernel $w(r) = \exp[-r^2/(2\sigma^2)]\sigma^{-2}(2\pi)^{-1}$ for which $q_k = 1/\sigma Q(L_k)$ with $Q(l) = \sqrt{-2\ln(1-l)}$ and compared the results with those obtained above with an exponential kernel. We did not observe significant differences with the exponential kernel, except at very high frequency. Such differences at high frequency are better seen with long-range correlated fields, as visible in Fig. 8: The Gaussian kernel leads to an underestimation of the spectrum at large wave numbers. This can be ascribed to finite numerical precision. With the Gaussian kernel, the eigenvalues $\hat{L}_k$ are distributed closer to 1, which is a singularity of $Q$. For this reason, a higher precision on $\hat{L}_k$ would be required in the Gaussian case.

### E. CFT illustrated with coarsening foam

To illustrate the cellular Fourier transform, we analyzed data for a two-dimensional foam shown in Fig. 1. We first obtained cell contours using the MATLAB imfill function which is based on morphological reconstruction. We ascribed pixels to cells; we calculated the discrete Laplacian from (13) by summing over pixels and then computed the harmonics. Figure 9 shows the spectrum of areal growth and cell density.

<figure><img src="image_10.png" /><figcaption>FIG. 7. Effect of kernel decay length. Spectra of stationary waves (a) and root-mean-square spectra of long-range correlated noise (b). The stationary waves are those represented in Fig. 2(a) and one realization of the random signal was shown in Fig. 2(g). The different colors correspond to different kernel decay lengths $\sigma$ of the exponential kernel. Black lines represent the Fourier transform of the input signals $f$ (theoretical curves).</figcaption></figure>

The spectrum for growth seems to be random and overall decay with wave number. Its resemblance with the spectrum shown in Fig. 6(c) suggests long-range correlations for areal growth in foams. Knowing that, in a coarsening foam, relative areal growth of a cell of area $S$ and number of neighbors $n$ is proportional to $(n-6)/(Se)$ with $e$ average thickness of films, this could be associated with long-range correlations in any of these parameters. The spectrum for cell density $1/S$ is rather flat, with a slight increasing trend at high wave numbers. This trend would be expected based on a consequence of the Aboav-Weaire and Lewis laws, according to which relatively small cells are in contact with relatively large cells (and conversely) leading to a short-range anticorrelation in cell size [43]. Although this analysis is preliminary and should be taken with caution, it indicates that CFT may help reveal new features of geometrically disordered materials.

## V. CONCLUSION

We considered here the harmonic decomposition of fields defined over disordered cellular media. We found that the classical fast Fourier transform and graph Fourier transform
<figure><img src="image_11.png" /><figcaption>FIG. 8. Effect of kernel type on the spectra. Log-log plot showing the mean-square spectra of the representations of fields with long-range correlations, as obtained with Gaussian (red) and exponential (blue) kernels ($\sigma = l_c$). Spectra were averaged over 1000 realizations. The black line represents the same root-mean-square spectra but deriving from the Fourier transform of the input signals $f$ (theoretical curve).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 9. Coarsening foam. Spectrum of areal growth in (a) and cell density in (b) for the foam shown in Fig. 1. Kernel: exponential, $\sigma = 7l_c$.</figcaption></figure>

have some shortcomings in this context. In brief, the FFT is defined on a square grid and does not account for disordered topology (variable number of neighbors), while the GFT does not account for disordered cell shape. We therefore built a more adequate harmonic decomposition, which we called cellular Fourier transform. It is based on the definition of a coarse Laplace operator and the use of an appropriate localized kernel. We found that the resulting harmonics generally behave as expected for an exponential kernel of decay length that is intermediate between cell size and domain size, though the CFT does not yield perfect peaks because of finite cell size. Another limitation of the CFT is that it is not well suited to the identification of wave vectors because it was designed for domains of irregular shapes. Overall, this harmonic decomposition is suited to disordered media that are divided into cells with variable sizes and irregular arrangements, such as foams, emulsions, granular materials, or biological tissues. As the definition of the harmonics does not depend on the coordinate system, our approach would also be applicable to non-Euclidean geometries, such as curved surfaces embedded in 3D and, in particular, biological thin tissues with complex 3D shapes.

Our method could be broadly useful for disordered media, even in the absence of subdivisions into cells. In some experimental situations like in fluid dynamics, it is possible to track landmarks such as particles to quantify their displacement. To a certain extent, and similarly to the Helmholtz decomposition of a vector field as the sum of a curl-free field and a divergence-free one, an equivalent point of view would be to consider a triangulation of the landmark distribution and to study the deformation and the rotation of the triangles. One could define invariants which do not depend on the translation relative to the reference frame or the coordinate system and directly apply the CFT to those invariants.

Finally, we believe that our approach can be used for cellularized media in all contexts where Fourier transforms are used. This includes statistical estimations (estimating fluctuations at different scales and their correlations), constructing a wavelet decomposition, or pattern recognition.

## ACKNOWLEDGMENTS

We thank Pierre Borgnat and Chun-Biu Li for fruitful discussions and Jérôme Duplat for sharing raw data of foam coarsening and for useful comments. This work was supported by the Université de Lyon through the program “Investissements d’Avenir” (ANR-11-IDEX-0007), and by the French National Research Agency through a European ERA-NET Coordinating Action in Plant Sciences (ERA-CAPS) grant (Grant No. ANR-17-CAPS-0002-01).

## APPENDIX: TO COMPUTE THE CELLULAR FOURIER TRANSFORMS

We here explain how to implement the cellular Fourier transforms of cellularized signals. This implementation comprises three steps: computing the Laplace operator, defining the Fourier harmonics, and computing the cellular Fourier transform of the cellular representation of a signal. In complement, we provide a MATLAB code (see Supplemental
Material [44]) that can be used to compute cellular Fourier transforms of a signal on a Voronoi tessellation.

1. Calculation of the discrete Laplace operator

We assume the tessellation to contain N cells, $\omega_i$, indexed by $i$ spanning $\{0, 1, \dots, N-1\}$. The geometry of a cell $\omega_i$, may be defined from an image, by the set of pixels that belong to $\omega_i$, or, alternatively, from a triangle mesh by the set of triangles that belong to $\omega_i$. A cell may also be defined by its contour, for instance by the set of pixels that belong to the contour or by the coordinates of the vertices of a polyhedron that defines the contour. It is then possible to compute the size $S_i$ of a cell (e.g., its area in 2D), typical cell scale $l_c = \sqrt{1/N \sum_{i=0}^{N-1} S_i}$, and the Laplace operator, whose elements are

$$ \bar{L}_{ij} = \delta_{ij} - \bar{W}_{ij}, \quad \text{with} \\ \bar{W}_{ij} = \sqrt{\frac{S_i \int_{\omega_i} \int_{\omega_j} w(|x-y|) d\mu(x) d\mu(y)}{S_j \int_{\omega_i} \int_{\Omega} w(|z-t|) d\mu(z) d\mu(t)}}. \quad (A1) $$

Here $w$ is the kernel function. In the paper, we generally used $w(x) = e^{-x/\sigma}$ with $\sigma = 7l_c$. The implementation of the integrals in the expression of $\bar{W}$ differs according to how the information on cellular geometry is stored. Any integral may be approximated by a sum over the elements (pixels, triangles) of the cell or, using the Stokes theorem to get a contour integral, by a sum over the elements of the contour. In the code that we provide, the functions `cft_voro` and `cft_laplace` generate random Voronoi tessellations and calculate their Laplace operator.

2. Calculation of the Fourier harmonics

The Fourier harmonics $e_k = \sum_{i=0}^{N-1} U_{ki} \Psi_i$ are computed by determining the transformation matrix $U$ from the basis

$\{\Psi_i\}_{i=0...N-1}$ of the representation space,

$$ \psi_i : \begin{array}{c} \Omega \to \mathbb{R}, \\ x \to \psi_i(x) = \begin{cases} 1/\sqrt{S_i} & \text{if } x \in \omega_i \\ 0 & \text{else.} \end{cases} \end{array} \quad (A2) $$

$U$ is defined as the right singular vectors of the Laplace operator. A singular value decomposition algorithm applied to $\bar{L}$ yields the left singular vectors $V$, the right singular vectors $U$, and the singular values $\hat{L}_k$:

$$ \bar{L}_{ij} = \sum_{k=0}^{N-1} \hat{L}_k V_{ki} U_{kj}, \quad \text{with} \quad \hat{L}_k = \hat{L}(q_k)/\hat{L}(0). \quad (A3) $$

To find the wave numbers $\{q_k\}_{k=0...N-1}$ associated with the harmonics, we use the singular value and the relation $\hat{L}_k = \hat{L}(q_k)/\hat{L}(0)$ where $\hat{L}$ depends on the functional form of the kernel $w$. For an exponential kernel $w(x) = e^{-x/\sigma}$, one has $q_k = 1/\sigma Q(\hat{L}_k)$, with $Q(l) = \sqrt{(1-l)^{-2/3} - 1}$. The harmonics are indexed so that their index grows with the wave number. The function `cft_harmonics` of our code computes the Fourier harmonics.

3. Calculation of the CFT of a signal

To calculate the cellular Fourier transform of a signal $f$, we compute its right coordinates in the basis $\{\Psi_i\}_{i=0...N-1}$ of the representation space, $\Phi_i = \langle \Psi_i \cdot f \rangle$. This done, the cellular Fourier transform is simply $\hat{\Phi}_k = \sum_{i=0}^{N-1} U_{ki} \Phi_i$, and the Fourier spectrum is obtained from $\hat{\Phi}_k$ and the wave number $q_k$. The functions `cft_cellularize_signal` and `cft_fourier` of our code compute signals in the basis $\{\Psi_i\}_{i=0...N-1}$ of the representation space and compute the Fourier transform.

[1] J.-F. Sadoc and N. Rivier, *Foams and Emulsions* (Springer Science & Business Media, New York, 1999), Vol. 354.

[2] D. L. Weaire and S. Hutzler, *The Physics of Foams* (Oxford University Press, New York, 2001).

[3] D. Weaire and M. A. Fortes, Stress and strain in liquid and solid foams, *Adv. Phys.* **43**, 685 (1994).

[4] L. J. Gibson and M. F. Ashby, *Cellular Solids: Structure and Properties* (Cambridge University Press, Cambridge, UK, 1999).

[5] S. Slotterback, M. Toiya, L. Goff, J. F. Douglas, and W. Losert, Correlation Between Particle Motion and Voronoi-Cell-Shape Fluctuations During the Compaction of Granular Matter, *Phys. Rev. Lett.* **101**, 258001 (2008).

[6] X. Trepat and E. Sahai, Mesoscale physical principles of collective cell organization, *Nat. Phys.* **14**, 671 (2018).

[7] E. K. Paluch, Biophysics across time and space, *Nat. Phys.* **14**, 646 (2018).

[8] I. Fortuna, G. L. Thomas, R. M. C. de Almeida, and François Graner, Growth Laws and Self-Similar Growth Regimes of Coarsening Two-Dimensional Foams: Transition from Dry to Wet Limits, *Phys. Rev. Lett.* **108**, 248301 (2012).

[9] S. Hilgenfeldt, S. A. Koehler, and H. A. Stone, Dynamics of Coarsening Foams: Accelerated and Self-Limiting Drainage, *Phys. Rev. Lett.* **86**, 4704 (2001).

[10] P. Mayer, H. Bissig, L. Berthier, L. Cipelletti, J. P. Garrahan, P. Sollich, and V. Trappe, Heterogeneous Dynamics of Coarsening Systems, *Phys. Rev. Lett.* **93**, 115701 (2004).

[11] J. Lambert, R. Mokso, I. Cantat, P. Cloetens, J. A. Glazier, F. Graner, and R. Delannay, Coarsening Foams Robustly Reach a Self-Similar Growth Regime, *Phys. Rev. Lett.* **104**, 248304 (2010).

[12] Y. Furuta, N. Oikawa, and R. Kurita, Close relationship between a dry-wet transition and a bubble rearrangement in two-dimensional foam, *Sci. Rep.* **6**, 37506 (2016).

[13] J. Duplat, B. Bossa, and E. Villermaux, On two-dimensional foam ageing, *J. Fluid Mech.* **673**, 147 (2011).

[14] A. Donev, F. H. Stillinger, and S. Torquato, Unexpected Density Fluctuations in Jammed Disordered Sphere Packings, *Phys. Rev. Lett.* **95**, 090604 (2005).

[15] A. J. Liu and S. R. Nagel, *Jamming and Rheology: Constrained Dynamics on Microscopic and Macroscopic Scales* (CRC Press, Boca Raton, FL, 2001).
[16] T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Jamming Transition in Granular Systems, Phys. Rev. Lett. 98, 058001 (2007).

[17] A. R. Abate and D. J. Durian, Approach to jamming in an air-fluidized granular bed, Phys. Rev. E 74, 031308 (2006).

[18] A. J. Liu and S. R. Nagel, Jamming is not just cool any more, Nature (London) 396, 21 (1998).

[19] G. Katgert, B. P. Tighe, and M. van Hecke, The jamming perspective on wet foams, Soft Matter 9, 9739 (2013).

[20] Y. J. M. Brechet and D. Weygand, Cellular structures in metallurgy, Foams and Emulsions (Springer, New York, 1999), pp. 151–174.

[21] T. Lecuit and P.-F. Lenne, Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis, Nat. Rev. Mol. Cell Biol. 8, 633 (2007).

[22] C. P. Gonatas, J. S. Leigh, A. G. Yodh, J. A. Glazier, and B. Prause, Magnetic Resonance Images of Coarsening Inside a Foam, Phys. Rev. Lett. 75, 573 (1995).

[23] A. J. Meagher, M. Mukherjee, D. Weaire, S. Hutzler, J. Banhart, and F. Garcia-Moreno, Analysis of the internal structure of monodisperse liquid foams by x-ray tomography, Soft Matter 7, 9881 (2011).

[24] W. Denk and H. Horstmann, Serial block-face scanning electron microscopy to reconstruct three-dimensional tissue nanostructure, PLoS Biol. 2, e329 (2004).

[25] J. B. T. M. Roerdink and A. Meijster, The watershed transform: Definitions, algorithms and parallelization strategies, Fundam. Inform. 41, 187 (2000).

[26] J. Lambert, I. Cantat, R. Delannay, R. Mokso, P. Cloetens, J. A. Glazier, and F. Graner, Experimental Growth Law for Bubbles in a Moderately “Wet” 3D Liquid Foam, Phys. Rev. Lett. 99, 058304 (2007).

[27] M. Saadatfar, A. P. Sheppard, T. J. Senden, and A. J. Kabla, Mapping forces in a 3D elastic assembly of grains, J. Mech. Phys. Solids 60, 55 (2012).

[28] P. B. de Reuille, A.-L. Routier-Kierzkowska, D. Kierzkowski, G. W. Bassel, T. Schüpbach, G. Tauriello, N. Bajpai, S. Strauss, A. Weber, A. Kiss et al., Morphographx: A platform for quantifying morphogenesis in 4D, Elife 4, e05864 (2015).

[29] A. Fruleux and A. Boudaoud, Modulation of tissue growth heterogeneity by responses to mechanical stress, Proc. Natl. Acad. Sci. USA 116, 1940 (2019).

[30] B. Helffer, Spectral Theory and its Applications (Cambridge University Press, Cambridge, UK, 2013), Vol. 139.

[31] L. Tan and J. Jiang, Digital Signal Processing: Fundamentals and Applications (Academic, New York, 2018).

[32] S. Bagchi and S. K. Mitra, The Nonuniform Discrete Fourier Transform and its Applications in Signal Processing (Springer Science & Business Media, New York, 2012), Vol. 463.

[33] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains, IEEE Signal Process. Mag. 30, 83 (2013).

[34] H. Zhang, O. van Kaick, and R. Dyer, Spectral methods for mesh processing and analysis, in Proceedings of Eurographics State-of-the-Art Report (Eurographics Association, 2007), Vol. 122.

[35] F. R. K. Chung and F. C. Graham, Spectral Graph Theory, CBMS Regional Conference Series in Mathematics (American Mathematical Society, Providence, RI, 1997), Vol. 92.

[36] D. A. Spielman, Spectral graph theory and its applications, in 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS '07) (IEEE, New York, 2007), pp. 29–38.

[37] M. Reuter, S. Biasotti, D. Giorgi, G. Patanè, and M. Spagnuolo, Discrete Laplace-Beltrami operators for shape analysis and segmentation, Comput. Graph. 33, 381 (2009).

[38] S. Rosenberg and R. Steven, The Laplacian on a Riemannian manifold: An introduction to analysis on manifolds (Cambridge University Press, Cambridge, UK, 1997), Vol. 31.

[39] M. Wardetzky, S. Mathur, F. Kälberer, and E. Grinspun, Discrete Laplace operators: No free lunch, in Symposium on Geometry Processing (Aire-la-Ville, Switzerland, 2007), pp. 33–37.

[40] H. A. Makse, S. Havlin, M. Schwartz, and H. Eugene Stanley, Method for generating long-range correlations for large systems, Phys. Rev. E 53, 5445 (1996).

[41] In the limit $\sigma \to 0$, the coarse Laplace operator converges toward $C\nabla^2 f$ in the bulk and toward $A(N \cdot \nabla)f$ at the boundaries where $N$ is the unit normal vector oriented inward $\Omega$ and the values of the constants are $C = -\int_{\mathbb{R}^n} w(|z|)z^2/2 d\mu(z)$ and $A = \int_{-\infty}^{+\infty} \int_0^{+\infty} w(\sqrt{x^2+y^2})y dx dy$. An eigenfunction $f$ of $\mathcal{L}$ associated to the eigenvalue $\lambda$ verifies therefore $C\nabla^2 f = \lambda$ in the bulk and $A(N \cdot \nabla)f = \lambda$ at the boundaries. Moreover in the same limit $\sigma \to 0$, the eigenvalue $\lambda$ is negligible at the boundaries, and the eigenfunctions of $\mathcal{L}$ are also the eigenfunctions of the Laplace operator with Neumann boundary conditions on $\Omega$.

[42] H. P. W. Gottlieb, Eigenvalues of the Laplacian with Neumann boundary conditions, ANZIAM J. 26, 293 (1985).

[43] S. N. Chiu, Aboav-Weaire's and Lewis' laws—A review, Mater. Charact. 34, 149 (1995).

[44] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.3.023036 for Matlab functions for the application of cellular Fourier transform to 2D Voronoi tessellations.