Cellular Fourier analysis for geometrically disordered materials

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 parameters and fields reflecting material properties and state are often assessed. In these media, it is difficult to quantify spatial variations of cell-scale fields, because space discretization based on standard coordinate systems is not commensurate with the natural discretization into geometrically disordered cells. Here we consider the spectral analysis of spatially varying fields. We built a method, which we call Cellular Fourier Transform (CFT), to analyze cell-scale fields, which includes both discrete fields defined only at cell level and continuous fields smoothed out from their sub-cell 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 biological tissue with complex geometry, opening the way to systematic multiscale analyses of material behavior.


Introduction
The past decades have seen a growing interest in geometrically disordered media [1] such as liquid and solid foams [2,3,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,9,10,11,12,13], fluctuations [14], jamming transition [15,16,17,18,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 3D. 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 2D and 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 are concerned with quantitative analyses of properties or fields defined on such segmented images.
An example of such medium is given in Figure 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, Figure 1) because growth requires landmarks (here, vertices) to be computed. These specificities 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 this difficulties, based on harmonic representation of signals defined on cellular media, which makes it possible to properly analyze the spectra of these signals. We term our approach Cellular Fourier Transform (CFT).
Spectral analysis decompose 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 irregular grids were based on FFT [32]. More recent approaches define 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 [35,36,33]. A graph may be endowed with an irregular geometry by ascribing a distance to each edge and -0.005 0.005 Areal growth rate in 2D foam Figure 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].
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. Some are related to the weak convergence of these operators towards the smooth Laplace-Beltrami operator in the limit of small mesh size [38] and the fact that discrete Laplacians on triangular meshes can not satisfy all desired natural properties [39]. Other reasons are related to the nature of the mathematical object that we consider: cellular tilings of the space (2D or 3D). In these tilings, cells may have complex shapes and 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) tilings of space. We start by presenting the geometry of the medium and how signals are represented. We define a coarse Laplace operator, applicable to signals with variations at sub-and supra-cell 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.

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 Ω of the n-dimensional space R n , divided into N nonempty subdomains {ω i , i = 0, 1, . . . , N − 1}, which we call cells. We consider a measure defined on R n , such as area in R 2 . Let µ and µ i be the measures (e.g. areas) of domain Ω and cells ω i , respectively: µ = N−1 i=0 µ i . We aim at analyzing a signal f defined over Ω (and assumed to have sufficient regularity for all the mathematical formulation to be well-posed). The smoothed version of the signal f is given by the piecewise constant function, where ω i dµ(y) (.) is the integral over the domain ω i and dµ(y) is the measure density (dµ(y) = dy for a Euclidean metric). We call f C the representation of f -it replaces f (x) by the average of f over the cell to which x belongs. We show in Figures 2 examples of signals f and their representations f C on a random Voronoï tessellation (for which the positions of the seeds have been chosen randomly, see Section 2.1). The representation of a function f C belongs to a vectorial space E of dimension N (the number of cells ω i ), which we call the representation space. A basis of E is the set of functions {ψ i } that vanish outside ω i and are defined by Given the usual scalar product f · g = Ω dµ(x) f (x)g(x) of two functions f and g, the basis of functions {ψ i } is orthonormal, The representation f C of f is also its orthogonal projection on the representation space E, In the following, we aim to define another orthogonal basis for the representation space 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 = i U ki ψ i , where U ki = e k · ψ i are the coefficients of the transformation matrix U between the two bases. Finding the harmonics of the representation space is equivalent to determine the unitary matrix U.

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 one of such 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 .

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 W which, to each function f defined on R n , associates The kernel w is an integrable real function and |x − z| is the Euclidean distance between points x and z of R n . We assume that w(r) has a maximum at r = 0 and vanishes when r → ∞, with a characteristic decay length σ. Like the discrete Laplacian on a grid, it averages the difference between the local field f (x) and its value on the neighborhood of x. In the limit where the lengthscale σ vanishes, W f ≃ C∇ 2 f , where ∇ 2 is the classical Laplace operator and the constant C = − R n w(|z|) z 2 /2 dµ(z). The operator W can be seen as a coarse version of ∇ 2 . In all generality, plane waves of wavenumber q, u q : x → exp(Iq · x), (with I 2 = −1) are eigenfunctions of W: . We will later consider the case whenL is positive and monotonously increasing on [0, +∞[ so that the eigenfunctions of W associated to a given eigenvalue are linear combinations of plane waves with all wavenumber having the norm |q|.

In bounded domain
We now consider a compact bounded domain Ω of R n and functions f defined on Ω, with L is a compact and self-adjoint operator; according to the spectral theorem, the eigenfunctions of L form an orthogonal basis of the associated Hilbert space [30]. If w is well localized around 0 (equivalently if its decay length σ is small with respect to size of domain Ω), then the integral in (5) can be approximated by an integral over the whole Euclidean space R n , as long as the position x is not too close to the boundariy of Ω. Therefore, Equation (4) implies that the eigenvalues of L are approximatelyL(|q|). Associated eigenfunctions are locally approximated by linear combinations of plane waves with wavenumbers of norm |q|, except close to the boundary; eigenfunctions may show boundary layers of width σ (see below). Note that, in order to reduced this boundary effect for eigenfunctions associated to the lowest eigenvalues, we defined L[ f ](x) in Eq. (5) using the difference f (x) − f (y), which enforces constant functions to be eigenfunctions of L and to be associated to the eigenvalue 0.

Harmonics of the representation space
In this section we discretize the problem and consider the representation of the coarse Laplace operator defined above. We introduce in section 1.3.1 the harmonics of the representation space and investigate in section 1.3.2 how they are related to the eigenfunctions of L and how their wavenumber can be computed. In section 1.3.3, we introduce a correction to reduce boundary effects.

Discrete Laplace operator
In order to build the representation of the coarse Laplace operator, we first consider the representation (the discretized version) w C (x, y) of w(|x − y|). It takes the form w C (x, y) = N−1 i=0 N−1 j=0 W i j ψ i (x)ψ j (y), {ψ i } being the basis of the representation space. The elements of the weight matrix W are given by The integral operator L C associated to the kernel For a function f C in the representation space, we may write We call L the discrete Laplace operator. It shows similarities with classical discrete Laplace operators [39], such as being symmetric, having positive weights W i j , and being positive semi-definite. 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 we consider here cellularized spaces on which functions are piecewise constant. The discrete Laplace operator L can be diagonalized and we propose to define the unitary matrix U introduced in section 1.1 from the eigenvectors of L, yielding where {L k , k = 0, 1, . . . , N −1} are the eigenvalues of L. Note that the eigenvector associated to eigenvalue 0 is ( instead of being (1, 1, . . . , 1) as for classical discrete Laplace operators. The harmonics {e k , k = 0, 1, . . . , N − 1} associated to U are the eigenfunctions of L C :

Eigenvalues
We consider an eigenfunction f of L: . Based on the definitions above, it is also such that its representation f C verifies, In the integral in (10), f (y) − f C (y) varies quickly, i.e. at the typical cell scale l c = (µ/N) 1/n . If the decay length σ of the kernel w is greater than or comparable to l c , then the integral in (10) is negligible and f C is, within a good approximation, an eigenfunction of L C . Following Eq. (9), the harmonics e k are the approximate representations of the eigenvectors associated to the eigenvalueŝ L k . Recall that we concluded in section 1.2 that each eigenfunction of L associated to eigenvalueL(|q|) is locally well approximated by a linear combination of plane waves with wavenumbers having the norm |q|. Altogether, the e k are appropriate harmonics of the representation space in the sense that they are locally approximated by a combination of plane waves of the same wavenumber q k given byL

Correcting boundary effects
The results obtained in Sections 1.
This rescaling only marginally changes the matrix L, though it breaks its symmetry. For this reason U is no longer defined as the eigenvectors ofL but as its right-singular vectors. This does not significantly affect U, except fo columns corresponding to cells at the edges of Ω. The singular valuesL k of the rescaled Laplace operator are therefore still related to wavenumbers via Eq. (11). If V are the left-singular vectors, we may write the components ofL as, From this analysis, we can deduce the spatial spectrum of a signal f . The k-th spectral coefficient is given bŷ and is associated to the wavenumber q k . We call this spectrum the Cellular Fourier Transform: It is the appropriate equivalent of the classical Fourier spectrum in R n for signals defined at cell level. 2 Implementation and results

Implementation
To test the Cellular Fourier Transform (CFT), we generated a domain and its partition into cells by generating the Voronoï 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 Fig. 3. Ω is not a perfect square but has a polygonal shape with a large number of edges and the number of cell neighbors is broadly distributed around 6. We tested several kernels w and characteristic decay lengths σ (see below). We numerically computed the integrals defining the kernel representation in Eq. (12) using Gauss quadrature. For tractability of numerical calculations, we used the approximation w(r > 6σ) = 0. We obtained the harmonics and associated wavenumbers using the singular value decomposition algorithm of MATLAB (Mathworks). We applied the spectral decomposition to a few deterministic and stochastic scalar functions defined over Ω. We present in Section 2.2 the results of our approach before testing in Section 2.3 how the parameters of the kernel influence the analysis. As an illustration, we analyze in Section 2.4 the foam coarsening data shown in Figure 1.

CFT applied to artificial fields
In this section, we use an exponential kernel, w(x) = σ −d exp(−r/σ). Following (11), the corresponding wavenumbers q k take values q k = 1/σ Q(L k ), where Q(l) = (1 − l) −2/3 − 1 is the inverse function ofL introduced in (4). Note that the wavenumber q k associated to the harmonics e k is not given by the square root ofL k , as would be for the Fourier transform in infinite space (though q k ∼ L k for smallL k ). These wavenumbers are shown in Figure 4. Figure 3 shows the first harmonics e k in Ω. These harmonics resemble the eigenfunctions of the Laplace operator with Neuman Indeed, the coarse Laplace operator converges to the classical Laplace operator with these specific boundary conditions in the limit σ → 0. The first harmonic, e 0 is for example constant as u 00 is. The following harmonics e 1 and e 2 correspond respectively to u 10 and u 01 respectively but their orientation deviates slightly toward diagonals. More generally, a given harmonic e k corresponds to a linear combination of eigenfunctions u mn for which π/a √ m 2 + n 2 ∼ q k . This linear combination is such that higher harmonics e k tend to have the same spatial periodicity in all directions.
We show in Figure 2 different fields (panels A-D) and their representations (E-H). We start with stationary waves (Fig. 2.A and E) of wavenumbers 2 √ 2πn/a with n = 1, 3, 5 and 10. We then consider a step function f (x, y) = −1 if y < 0 and f (x, y) = 1 if y ≥ 0 ( Fig. 2B and F). Then, we built a long-range correlated random field ( Fig. 2C and G) using the Fourier filtering method detailed in [41]; it gave the image Fig. 2C for which the 2 points correlation function for which pixel intensity I i I j decays with the distance d i j between pixels like I i I j ∝ 1/d i j . Finally, we considered a white noise (Fig. 2D and H).
In Figure 5 we show the spectra of the fields plotted in Figure 2. The spectra of the stationary waves representations are indeed maximal close to the wavenumbers of the initial waves ( (Fig. 5bf A). Also as expected, the spectrum of the step function representation is peaked at 0 (Fig. 5bf B). The spectra of field with long-range correlations (Fig. 5C) and of the white noise (Fig. 5D) 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 both white noise and long-range correlated field as above. We estimated the average spectrum, which should give an estimate of the Fourier transform of the correlation function, and the spectra of representations behave as expected: constant for white noise and power-law decay for field with long-range correlations (Fig. 5E).

Sensitivity of the CFT to the kernel
We first tested the influence of the kernel decay length, σ, of the exponential kernel. We show in Figure 6 spectra obtained with different values σ. We see that when σ is small, spectra are shifted towards higher wavenumbers. Optimizing the kernel requires taking into account domain size and numerical precision at which the singular value decomposition is performed. The cell size l c should be a lower bond for σ while the domain size and the numerical precision prevent σ to be too large. Optimal value for σ 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 wavenumbers corresponding to maxima of spectra in 6A. and the wavenumber of the corresponding stationary waves. We found σ ≃ 7l c , which is intermediate between cell size and domain size.
We calculated the harmonics of the same domain Ω using a Gaussian kernel w(r) = exp(−r 2 /(2σ 2 ))σ −2 (2π) −1 for which q k = 1/σQ(L k ) with Q(L k ) = −2 ln(1 −L k ) 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 Figure 7: The Gaussian kernel leads to an underestimation of the spectrum at large wavenumbers. This can be ascribed to finite numerical precision. With the Gaussian kernel,L k are distributed closer to 1 which is a singularity of Q. For this reason, a higher precision onL k would be required in the Gaussian case.  Figure 5: CFT of the cellularized fields shown in Figure 2. A-D Spectra of the cellularized field shown in Figure 2E-H, respectively. E Root mean square of the spectra of white noise (red) and of random fields correlated at long range (blue); the average is over 1000 realizations. Dashed lines in A are at the wavenumber of the initial fields. The dark lines in C, D and E represent the theoretical curves given by the Fourier transforms of the continuous signals f (theoretical curve). Kernel: exponential, σ = 7 l c . Stationary waves Random signals Figure 6: Effect of kernel decay length. Spectra of stationary waves A and root mean square spectra of long-range correlated random signals B. The stationary waves are those represented in Fig. 2A and one realization of the random signal was shown in Fig. 2G. The different colors correspond to different kernel decay length σ of the exponential kernel. Black lines represent the same root mean square spectra but deriving from the Fourier transform of the input signals f (theoretical curve). Gaussian v.s. ex Figure 7: 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 (σ = 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).

CFT illustrated with coarsening foam
To illustrate the Cellular Fourier Transform, we analyzed data shown in Figure 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 (12) by summing over pixels and then computed the harmonics. Figure 8 shows the spectrum of areal growth, which seems to be random and overall decay with wavenumber. Its resemblance with the spectrum shown in Figure 5I suggests long-range correlations for areal growth in foams, which might reflect their elastic behavior at short timescales. Although this analysis is preliminary, it indicates that CFT may reveal new features of geometrically disordered materials.

Conclusion
We built here a harmonic decomposition for signals defined over cellularized media, which we call Cellular Fourier Transform. It is based on the definition of a coarse Laplace operator and the use of an appropriate kernel. We found that the resulting harmonics behave as expected for an exponential kernel of decay length that is intermediate between cell size and domain size. 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 the translation relative to the reference frame or the coordinate system and directly apply the analysis 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. by the French National Research Agency through a European ERA-NET Coordinating Action in Plant Sciences (ERA-CAPS) grant (ANR-17-CAPS-0002-01).