Minkowski Tensors of Anisotropic Spatial Structure

This article describes the theoretical foundation of and explicit algorithms for a novel approach to morphology and anisotropy analysis of complex spatial structure using tensor-valued Minkowski functionals, the so-called Minkowski tensors. Minkowski tensors are generalisations of the well-known scalar Minkowski functionals and are explicitly sensitive to anisotropic aspects of morphology, relevant for example for elastic moduli or permeability of microstructured materials. Here we derive explicit linear-time algorithms to compute these tensorial measures for three-dimensional shapes. These apply to representations of any object that can be represented by a triangulation of its bounding surface; their application is illustrated for the polyhedral Voronoi cellular complexes of jammed sphere configurations, and for triangulations of a biopolymer fibre network obtained by confocal microscopy. The article further bridges the substantial notational and conceptual gap between the different but equivalent approaches to scalar or tensorial Minkowski functionals in mathematics and in physics, hence making the mathematical measure theoretic method more readily accessible for future application in the physical sciences.

The morphology of complex spatial microstructures is often classified qualitatively into types such as cellular, porous, network-like, fibrous, percolating, periodic, lamellar, hexagonal, disordered, fractal, etc. Various quantitative measures of morphology have been defined often applicable to one specific type only, for example moments of the distributions of angles of tangent vectors with a fixed specified direction as anisotropy characterization of a network structure. Apart from the concept of correlation functions, few measures are defined sensibly and robustly for all types. In this article, we describe the class of Minkowski tensors (MT) that apply generically to any type of structure, which contains two or more phases separated by a well-defined interface (e.g. porous media, foams, trabecular bone, granular material in Fig. 1). The MT are defined as integrals of powers of normal and position vectors and surface curvatures (or curvature measures). Because of their tensorial nature they are explicitly sensitive to anisotropic and orientational aspects of spatial structure. Robust indices of intrinsic anisotropy and alignment can be derived. Figure 1 shows examples of systems where subtle anisotropy of the spatial structure influences the physical properties and to which the analysis of this article is applicable. The scope of this article is to supply a thorough description of an algorithm to compute MT of bi-phasic structures, when the phase-separating interface is given by a triangulation. Two examples illustrate the Minkowski tensor analysis. For further applications of MT for pattern analysis we refer to refs. [1]- [5], [5]- [7]. Scalar Minkowski functionals (MF) are well established as succinct descriptors of morphology and spatial structure for various physical processes [10]. These integral geometric measures have been applied to disordered porous materials [11] and are relevant to flow phenomena therein [12], to nanoscale microstructures in copolymers [13], to the dewetting dynamics of thin films [14], and to Turing patterns [15]. They have also been shown to be the most pertinent morphological quantities on which the thermodynamic properties of simple fluids near curved solid interfaces depend [16], [17]. The mathematical theory of MF and their generalizations has been comprehensively developed in the context of integral geometry [18]- [21], with several aspects shared also with the discipline of mathematical morphology [22]- [24].
MF as scalar quantities are not sensitive to features of the morphology which relate to orientation or directional Examples of systems with anisotropic spatial structure. (a) A microphase-separated copolymer film aligns under the influence of an external electric field (adapted from [8], reproduced by permission of The Royal Society of Chemistry). (b) Closed-cell metal foam (image courtesy of M. Saadatfar [9]). (c) Structure of trabecular bone (image courtesy Alan Boyde). (d) Packing of ellipsoids as systems for anisotropic granular matter. anisotropy, since motion-invariance is one of their defining properties. Therefore they do not provide quantification of anisotropy that is relevant to study the direction-dependence of physical processes, such as elastic properties or permeability of anisotropic porous or microstructured materials or systems with external fields. This motivates their generalization to tensorial quantities. MT have already been shown to be the relevant morphological descriptors for a density functional theory of fluids of non-spherical particles [25] and of DNA conformations [26], and of a simple model for transport with molecular motors [27]. They have also been used, in 2D, as morphology descriptors of arrangements of neuronal cells [28], galaxies [29], and Turing patterns [6]. The mathematical discipline of integral geometry has proven statements regarding continuity and completeness similar to the Hadwiger theorem in the scalar case [30]- [33]. However, an algorithm for the computation of the MT applicable to experimental 3D data -a prerequisite for their use as shape indices for experimental data -has thus far been lacking.
A primary application of rank-2 MT is the quantitative analysis of the degree of intrinsic anisotropy of materials with complex spatial structure. Scalar measures of anisotropy are easily derived as eigenvalue ratios of the MT. Alternative methods for the characterization of anisotropy and alignment exist. Fourier transforms are a common way to characterize anisotropy, and have been applied e.g. for trabecular bone [34], for electrodeposited patterns [35], for fiber systems [36], and for structured polyethylene mats [37]. Related methods based on correlation functions are also known [38], [39]. Fourier methods that analyze the amplitude of the Fourier transform of a gray-scale image in polar coordinates can also quantify alignment, e.g. of copolymer films in electric fields [8]. Anisotropy indices derived from the normal vector distribution of a given shape, similar to the MT, have been used to describe the shape anisotropy of simulated 3D foam cells [40], [41] and liquid interfaces [3], [42]. An anisotropy measure applicable to porous media is derived from the directional variations of average chord lengths. For a binary composite, i.e. consisting of a solid and a void phase, a chord is a segment of an infinite straight line that is fully contained in one of the two phases. Analyses of chord lengths and the derived mean intercept length ellipsoid are used for the investigation of the microstructure of bone [43]- [48], see also ref. [49] for a comparison of anisotropy measures based on mean-intercept length, star-volume and star-length distributions. Deformations of cellular or granular material have recently been quantified using the so-called texture tensor C, defined as the sum C ij := l i l j over a subset of link vectors l in the structure [50], [51]. The texture tensor can be used to characterize anisotropy, e.g. for Antarctic ice crystals [52] and liquid foam cells [53]. Further anisotropy measures are based on the Steiner compact [54], wavelet analysis [55], the orientation of volumes [56] or star-volumes [57]. Twodimensional equivalents of the anisotropy measures discussed in this article have previously been used for the analysis of the shape of neuronal cells [28] and galaxies [29], and are discussed in detail in [6]. (Parts of the analysis of this article (a) A body K with bounding surface ∂K; (b) a convex polytope P and its parallel body (or dilation) P ⊎Bǫ; (c) a subset of a topologically more complex body based on Craig Marlow's painting "White Spirits" [59]. The latter demonstrates a common experimental situation, namely that the given body K represents a finite subset of a larger body K + , here the body K clipped to the window of observation T , K = K + ∩ T . K is assumed to be a representative subset of the larger body, which allows for the estimation of intrinsic shape features of K + . If only K, but not K + , is available for analysis particular care must be taken w.r.t. those bounding surface patches of K that result purely from taking the subset, i.e. those that are on the boundaries of the window of observation.
represent the thesis work in ref. [58]). The paper is organized as follows: Section I provides an overview of the theory of MF and MT, based on their definition by surface integrals which is widely used in the physical sciences. To bridge the gap in notation between the physics and the mathematics literature this section also includes a discussion of the alternative definition based on measure theory. Section II describes algorithms to compute MT for triangulated bodies. Section III describes anisotropy measures derived from the MT and illustrates their application to two experimental data sets. The appendix provides analytic expressions for the MT for some simple geometric shapes.

I. DEFINITION AND FUNDAMENTAL PROPERTIES
Scalar MF and MT can be used as shape measures to quantify the shape or form of an object. They can be defined in two largely equivalent ways. In the physical sciences, a definition based on curvature-weighted integrals of position or surface normal vectors over the object's bounding surface is commonly used, and forms the basis of the numerical algorithms derived in this article. An alternative, more fundamental definition is provided by the measure theoretic approach of integral geometry (see I-B).
The object, also referred to as body, whose shape can be characterized by MT or MF is denoted by K. Assuming that K is a compact set with nonempty interior embedded in Euclidean space Ê 3 and is bounded by a sufficiently smooth surface ∂K, we define MF of K as in space-dimension d = 3 and with ν = 1, 2, 3. The scalar functions G ν are G 1 = 1, the mean curvature G 2 = (k 1 + k 2 )/2 and the Gaussian curvature G 3 = k 1 ·k 2 of the bounding surface ∂K (k 1 , k 2 are the principal curvatures on ∂K as defined in differential geometry); dV is the infinitesimal volume and dA the scalar infinitesimal area element. This definition naturally applies to both convex and non-convex bodies of arbitrary topology with a sufficiently smooth bounding surface. The prefactor is chosen such that for a sphere B R with radius R the scalar MF are W ν (B R ) = κ 3 R 3−ν where κ 3 = 4π/3 is the volume of the 3-dimensional unit sphere 1 . The MF W 2 and W 3 (in 3D space) are not properly defined by eq. (1) for bodies with sharp edges or corners, due to singularities of the mean and Gaussian curvatures G 2 and G 3 . However, for convex bodies, consideration of a parallel or dilated body in the limit of vanishing thickness provides a robust definition, by use of the Steiner formula. The Steiner formula states that, for a given convex body K, the MF of the parallel or dilated body (K ⊎ B ǫ ) are a polynomial in ǫ with coefficients proportional to the MF of K; K ǫ := K ⊎ B ǫ := {x 1 + x 2 : x 1 ∈ K, x 2 ∈ B ǫ } is the parallel or dilated body of K (see Fig. 2). Specifically for the volume one finds Sharp edges and vertices of K correspond to cylindrical or spherical segments on the bounding surface ∂(K ⊎ B ǫ ) of the parallel or dilated body. The bounding surface is sufficiently smooth for ǫ > 0 and the limit ǫ → 0 of W ν (K ⊎ B ǫ ) converges to W ν (K). It is further necessary to define MF and MT for certain non-convex bodies (with or without positive reach, see e.g. [21, Note 1 to Section 5.3 and refs. therein]); this is achieved below by exploiting an additivity relationship. A further discussion for non-smooth bodies can be found in section II. MT are symmetric tensors (that is, invariant under index permutation), which are generated by symmetric tensor products of position vectors x and normal vectors n of ∂K. The dyadic (or tensor) product of two vectors a and b is (a⊗b) ij := a i b j . Let now a and b be symmetric tensors of rank r and s resp., then the symmetric tensor product is defined as where S r+s is the permutation group of r + s elements. We use for two tensors a and b the shorter notation a 2 := a⊙a = a ⊗ a and ab := a ⊙ b. For example, if a and b are both vectors the symmetric tensor product is ij of the rank 2 tensor ab. 1 Other normalizations are common in the literature too. The kinematic formulae [20], [21] have particularly simple coefficients if the normalization Mν (K) = κ d−ν Wν (K)/(κν κ d ) is used, with the volume of the ndimensional unit ball κn := π n/2 /Γ(n/2 + 1). In the mathematical literature, in d-dimensional Euclidean space the normalization V d−ν (K) = d ν Wν (K)/κν is frequently used, and the V d−ν are called the intrinsic volumes of K. In three dimensions, the set of MF thus consists of the volume W 0 = M 0 = V 3 , the surface area 3W 1 = 8M 1 = 2V 2 , the integrated mean curvature 3W 2 = 2π 2 M 2 = πV 1 , and the Euler characteristic The MT of rank two are defined as with ν = 1, 2, 3 and (r, s) = (2, 0), (1,1) or (0, 2). For ease of notation, we set W r,s 0 := 0 for s > 0 and W r,s ν := 0 if ν < 0 or ν > 3. For a three-dimensional body, this definition yields 10 MT (not counting the ones that vanish by definition for all bodies).
MT of rank one (called Minkowski vectors) are defined by W 1,0 0 := K x dV and by W 1,0 ν := 1 3 ∂K x dA for ν = 1, 2, 3. The prefactors are chosen such that, for a sphere centered at C, the so-called curvature centroids W 1,0 ν /W ν are equal to C. Note specifically that W 0,1 1 /W 0 is the center of mass (assuming a solidly filled body of constant density). Formally, vectors W 0,1 ν proportional to ∂K n dA for ν = 1, 2, 3 are also defined, however they vanish for any body (with a closed bounding surface) [32].
MT are isometry covariant, that is their behavior under translations and rotations is given by where K ⊎ t is the translation of K by the vector t, U K is a rotated copy of K, andÛ r+s denotes the corresponding rotation tensor for a rank-(r + s) tensor: In this expression, U ij is the conventional orthogonal 3 × 3 transformation matrix associated with the rotation U . For r = 0, equation (6) gives W r,s ν (K) = W r,s ν (K ⊎ t). A tensor that fulfills this relation for all K is called translation invariant. Genuinely translation covariant tensors fulfill (6) but not the translation invariance condition. For the sake of brevity, we will use the term translation covariant to denote specifically the genuinely translation covariant tensors. All are translation invariant due to the envelope theorems of Müller [32]. Translation invariance is important whenever a natural choice for the origin is not available.
All MF and MT are homogeneous, i.e. they fulfill the homogeneity relation W r,s ν (λK) = λ 3+r−ν W r,s ν (K).  without further assumptions (technically, such bodies do not represent sets of positive reach [21], [60]). An extension of the definition of MF and MT to finite unions of convex bodies is achieved by exploiting a property called additivity if K, K ′ are all convex. In general, the union (K ∪ K ′ ) of two arbitrary convex bodies K and K ′ is not convex while the intersection (K ∩ K ′ ) is convex. Since W r,s ν are continuous functionals, eq. (9) can be used (see Groemer's extension theorem [21]) to define the MF and MT for bodies that are unions of a finite number of convex bodies (such sets are called polyconvex).
MT of rank (r + s) with r + s > 1 are not completely linearly independent, i.e. they do not contain independent shape information [31], [61]. For rank-2 MT and d = 3 the linear dependences are valid for any polyconvex body K in Ê 3 and ν = 0, 1, 2, 3; in particular, it follows that QW 3 = 3W 0,2 3 . Specifically, for ν = 0 one obtains a special case of the Gauss' theorem with real coefficients ϕ r,s,p ν that do not depend on the convex body K. The coefficients ϕ r,s,p ν vanish unless r + s + 2p = 2. Starting from eq. (12) and using the linear dependences among the basic tensor valuations, ϕ can be expressed in terms of linearly independent basic tensor valuations which form a basis of the corresponding vector space. The vector space of continuous, isometry covariant tensor valuations of rank two in Ê 3 has dimension 10. A particular basis of this vector space consists of the six tensor valuations and W 0,2 2 , which contain important independent shape information, and of the four tensor valuations QW ν , ν = 0, . . . , 3. A summary is provided in Table I.
Since ϕ is continuous and additive on convex bodies, it can be extended as an additive functional to finite unions of convex bodies. For this additive extension, eq. (12) remains valid since the right-hand side is a linear combination of additive functionals. It should be emphasized, however, that although all these functionals are continuous on the space of convex bodies, they are not continuous on the space of finite unions of convex bodies, see the example in [6, Fig. 3].

A. MT of convex polytopes
It is instructive to illustrate the MT for convex polytopes P and to point out similarities to the tensor of inertia. We use here the letter P , rather than K, to denote a body whose bounding surface is a polytope, i.e. consisting of piece-wise linear facets. For a polytope, the tensors W 2,0 ν characterize the distribution of mass if the cell is a solid cell (W 2,0 0 ), a hollow cell (W 2,0 1 ), a wire frame (W 2,0 2 ) and a cell consisting of points at the vertices only (W 2,0 3 ); in the last two cases, however, this distribution of mass is weighted with certain exterior angles (see Fig. 3).
The tensor of inertia Á, defined by Á ij = K −x i x j + δ ij |x| 2 dV , is a measure of the mass distribution of a (homogeneous) body K, relevant for the relationship between a rotation and the resulting moment. As Á is not translation-covariant, it is not a linear combination of the MT. However the simple relationship Á(K) = −W 2,0 0 (K) + tr W 2,0 0 (K) Q holds for arbitrary K; as above Q is the unit tensor of rank two (metric tensor) of Ê 3 . This illustrates that W 2,0 0 (K) is also a measure of the distribution of mass if K is a homogeneously filled solid. Similarly, W 2,0 1 (K) characterizes the mass distribution if K is hollow and bounded by an infinitesimally thin surface sheet. This interpretation of W 2,0 0 and W 2,0 1 is valid for bodies K bounded by arbitrary surfaces (not just polytopes).
For a polytope P , the tensor W 2,0 2 (P ) reduces to a line integral over the edges of the polytope (as the mean curvature vanishes on the flat facets, see also section II) and is hence related to a mass distribution if P is given by a wire frame with wires along the edges. However, imposed by the requirement of additivity, the weight of the wire cannot be uniform but must be proportional to the mean curvature along the edge (i.e. the dihedral angle). Similarly, the tensor W 2,0 3 (P ) reduces to a sum of point contributions, as the Gaussian curvature G 3 of P vanishes except at the vertices of the given polytope P . Due to additivity again, these vertices are weighted with the Gaussian curvature G 3 .

B. Definition based on fundamental measure theory
This section describes the alternative (and in some sense more fundamental) definition of MF and MT based on integral geometry and fundamental measure theory. The purpose of this section is to bridge the gap between the mathematics and physics literature on MF and MT. However, its content is not required for the numerical approaches to MF and MT described in section II.
In integral geometry, the definition of MF and MT is based on so-called support measures (sometimes called generalized curvature measures) that can be thought of as local versions of the scalar MF [18], [20], [21], [62].
If support measures are available, then the MT for convex (or more general) bodies are obtained by integrating tensor functions with respect to these measures. Here we describe the approach for convex sets. The idea underlying the introduction of support measures for convex sets is to generalize the notion of a parallel set (or "dilation", see Fig. 4 Construction of a local parallel set. (a) Definition of the normal field over an arbitrary convex body K at its boundary ∂K. (b1-b2) Local parallel set Mǫ(P, η) (illustrated for a polytope P ) where only points are considered for which the normal direction is in a prescribed subset S of the unit sphere.
A definition of the local parallel set that also applies to convex bodies K without smooth boundary is given in the following: We define p K (x) as the unique point in K which is nearest to a given point x ∈ Ê d . This defines a continuous map p K (·) : Fig. 4). This construction achieves a definition of nearest points (reminiscent of the Euclidean distance map [63]) and surface normals that is also well-defined for points x whose nearest point on K is a sharp corner (where hence, the tangent plane is not well-defined and hence the conventional differential geometric definition of the surface normal does not apply). Now we shall derive a local Steiner formula and support measures, following ref. [21] to obtain local support measures. Then the definition of MF and MT follows directly as a special case [31], [61]. The intuitive idea underlying the definition of a local parallel set is described in two steps. First, we specify a region L ⊂ Ê d and some ǫ > 0. Then we consider all points x ∈ Ê d \ K which have distance d K (x) at most ǫ from K and for which p K (x) ∈ L. Second, if K has points (such as at sharp corners or edges) where the exterior surface unit normal vector is not unique, then it is natural to restrict the points in this local outer parallel set further by also requiring n K (x) to lie in a prescribed set S ⊂ Ë d−1 . As an example, consider the polytope P in Fig. 4 (b1); for the definition of the local parallel set (shaded dark) it is necessary to specify which angular fraction of the wedges should be part of the local parallel set. This motivates the definition of the local parallel set by specification of the spatial region L and by the subset sphere. This two step procedure can be merged and slightly extended to the following general definition.
For given ǫ > 0 and η ⊂ Ê d × Ë d−1 , the local parallel set of K defined by global outer parallel set (K ⊎ B ǫ ) \ K (see Fig. 4 The volume of this local parallel set of a convex body K is V d (M ǫ (K, η)). A fundamental result in integral geometry, known as the local Steiner formula [21], [62], states that the where κ n := π n/2 /Γ( n 2 +1) is the volume of an n-dimensional unit ball and Λ ν (K, η), ν = 0, . . . , d − 1, are certain real coefficients that depend on K and η, but not on ǫ.
For eq. (14) to be true for all ǫ > 0, it is crucial that K is convex [64]. Equation (14) is easily confirmed (and evaluated) for a convex polytope P . In this case, the set (K ⊎ B ǫ )\P can be decomposed in an elementary way into wedges over the faces F of P as indicated by Fig. 4 (b1). Let F ν (P ) denote the ν-dimensional faces of the polytope P 3 and the set of unit normal vectors assigned to F ∈ F ν (P ) and relint F the relative interior of F 4 ; see Fig. 5. The contributions to V d (M ǫ (K, η)) coming from each of these wedges where ω ν = ν κ ν is the surface measure of the (ν − 1)dimensional unit sphere and dH ν is the Hausdorff measure of dimension ν [67].
Since an arbitrary convex body K can be approximated by polytopes, eq. (14) can be derived by a continuity argument. In fact, Λ ν (K, η) can be expressed as a linear Therefore, the properties of the local parallel volume V d (M ǫ (K, η)) (in particular additivity, weak continuity) are also available for Λ ν (K, η) [62, p. 202]. In particular, K → Λ ν (K, η) is an additive functional for fixed η. That is, for convex bodies K and K ′ This property is called σ-additivity of Λ ν (K, ·). Weak continuity of Λ ν means that for every sequence of convex bodies K (µ) (µ ∈ AE), with K (µ) → K and every continuous function holds. Note that this does not imply are called support measures and are determined as coefficient measures of the Steiner formula (14). They are local versions of the classical MF W ν , since If K is sufficiently smooth, then where G ν (x) is the (ν − 1)-th (normalized) elementary symmetric function of the principal curvatures of ∂K at x. That is in three dimensions, these are 1, the mean curvature and Gaussian curvature, respectively. 1{·} is the characteristic function, which is evaluated to one if · is true and 0 otherwise. For general dimensions and ν ≤ d − 1, eq. (21) holds for all convex bodies K.
Having introduced the support measures as local versions of the scalar MF, it is easy to define the MT for a convex body K by hence we obtain Φ r,s ν (K) by integrating the tensorial function x r n s with respect to the measure Λ ν (K, ·) over N (K) ⊂ Ê d × Ë d−1 . If K is a polytope, for d = 3 this yields equation (48) (up to a different normalization). If K is smooth, we obtain eq. (5) (up to a different normalization and indexing scheme), i.e.
The notation Φ µ,r,s or Φ r,s µ for the MT in eq. (22) is preferred in the mathematical literature and differs from the notation W r,s ν in eqs. (4)(5) only by a different indexing scheme and a different normalization. In Ê 3 , i.e. for d = 3, the functionals Φ r,s µ and W r,s ν are related by for ν = 1, . . . , d, and The additivity and continuity properties of the support measures immediately yield the corresponding properties of the MT. This approach also shows that if it is possible to define support measures for a class of sets, then the corresponding tensor valuations can be defined by eq. (22).
Since the theory of support measures is well-developed [21], [62], the measure theoretic approach outlined above has some advantages over the differential geometric approach.
As a simple illustration, let us explain why W 1,1 d−ν is translation invariant for ν = 0, . . . , d − 1. Observe that by translation covariance of the support measures It is a basic property of the measures Λ ν (K, ·) that they are centered at the origin in the sense that n Λ ν (K, d(x, n)) = 0 [32], which yields the assertion.
A natural and useful extension that is suggested by general measure theory is to introduce local tensor valuations by restricting the integration on the right-hand side of eq. (22) to measurable subsets of Ê d × Ë d−1 .

II. BODIES BOUNDED BY TRIANGULATED SURFACES
We describe an exact algorithm for all independent scalar, vectorial and rank-2 MT of bodies bounded by piece-wise linear (i.e. triangulated) surfaces. Henceforth such bodies, convex or non-convex, are called polytopes P and their bounding surface is the triangulation F . Triangulations are commonly used as discrete approximations of smooth surfaces. The continuity of the MF and MT guarantee the convergence of the formula for the triangulations to the MF and MT of the smooth body.
The formulae are derived for convex bodies with triangulated bounding surfaces by considering parallel bodies P ǫ of convex polytopes P (that is P ǫ has a continuous normal fields and finite curvatures for ǫ > 0 and a well defined limit as ǫ → 0). By application of the additivity relation these formulae are then shown to be valid also for bodies that are not convex. As the key results of this article-explicit formulae for the computation of MT of convex and non-convex polytopesare summarized in table II.
Consider a polytope P in Ê 3 with piece-wise linear bounding surface ∂P ≡ F . Without loss of generality the linear facets may be assumed to be triangles. 5 The set of all triangular patches of ∂P is F 2 , the set of oriented edges is F 1 and the set of vertices is F 0 . We assume a doubly connected edge list (DCEL [68], also called half edge data structure [69]), that is, every edge which is shared between two triangles T and T ′ is a double-edge consisting of two oriented edges e (being part of T ) and e ′ (part of T ′ ), constituting an unambiguous assignment of each edge to a triangle. Each oriented edge is assigned to its previous edge e previous and its next edge e next . The remaining ambiguity in the edge orientation is lifted by requiring the triangle normals n T = (e previous × e)/ e previous × e to point out of the body P . Thus we can uniquely assign to each oriented edge e a triangle T with vertices v 1 , v 2 v 3 (see Fig. 6) and its normal vector n T .
The parallel body construction is illustrated by Fig. 2 (b). For an arbitrary body P the parallel body P ǫ with thickness ǫ > 0 is defined as P ǫ := P ⊎ B ǫ . For a convex polytope P (whose bounding surface has a discontinuous normal field) the bounding surface ∂P ǫ has a continuous normal field. The curvatures are patch-wise constant: G 2 = G 3 = 0 on the planar patches, G 2 = (2ǫ) −1 and G 3 = 0 on the cylindrical patches corresponding to polygon edges, and G 2 = 1/ǫ and G 3 = 1/ǫ 2 on the spherical patches corresponding to polytope vertices. For convex polytopes, the MT are defined as the surface integrals of eq. (5) evaluated on ∂P ǫ in the limit ǫ → 0. The result thus obtained is equivalent with eq. (22), see also eq. (48).
Second column: MF and MT for bodies with smooth boundary ∂K. The mean and Gaussian curvatures on ∂K are G 2 and G 3 , respectively. Third column: MF and MT for a triangulation. The set of facets of the triangulation F is F 2 , the set of oriented edges is F 1 (in DCEL structure, see text) and the set of vertices F 0 . The subset of triangles that contain the vertex v is denoted by F 2 (v). The nomenclature for triangulated surfaces is defined in Fig. 6. The vertices of an edge e or a triangle T are denoted v 1 , v 2 and v 3 , respectively. |T | is the area of T ∈ F 2 , C T its center point (v 1 + v 2 + v 3 )/3 and the tensors I T and J T are given in eqs. (33) and (35) and table III. Ce is the center point of edge e; Ce = (v 1 + v 2 )/2. i, j, k ∈ {x, y, z}. (see also [58]).  Fig. 7. Subdivision of a non-convex body into convex sub-bodies, P = P 1 ∪ P 2 ∪ P 3 ∪ P 4 ∪ P 5 . Note that for the computation of MT the segments P 2 and P 4 need to be taken into account, though their volume measure is 0: The calculation of the volume of a polytope P can be transformed into a surface integral by Gauss' law (see eq. (11)) [70]. With div x = div(x, y, z) t = 3 one obtains where dA = n dA denotes the oriented infinitesimal area element and ·, · the scalar product.

B. Surface area W 1 and integral mean curvature W 2
The surface integral is a sum over triangles and is easily evaluated yielding the formulae in Tab. II. This result is independent whether P is convex or not. The surface area W 1 (P ) of ∂P is simply the sum of triangle areas.
Expressing W 2 (P ) as the limit of vanishing parallel distance ǫ of W 2 (P ǫ ) of the parallel body P ǫ , W 2 (P ) = lim ǫ→0 W 2 (P ǫ ), the contributions of facets vanish because the mean curvature of a flat face is zero. The contribution of the spherical patches S corresponding to vertices vanishes because the integral over a spherical patch S can be parametrized in spherical coordinates by S 1 ǫ ǫ 2 sin θ dϕ dθ which vanish as ǫ → 0. The remaining contribution of the edges is located at the cylindrical patches of ∂P ǫ and is given in polar coordinates by where |e| is the length of edge e and α e the dihedral angle, i.e. the angle between the surface normals of the two facets adjacent to e. For a convex body, all edges have a dihedral angle 0 ≤ α e ≤ π; see also Fig. 6. Note that F 1 is the set of oriented edges, i.e. the edge shared by two triangles is represented by two distinct oriented edges, which explains the factor 1/2 in front of the integral.
Eq. (28) remains valid even if P is not convex, as is shown by exploiting additivity: A non-convex polytope P can be decomposed into a set of convex polytopes by cutting along the symmetric bisector planes of all concave edges (that is, −π < α e < 0), see Fig. 7. For a concave edge e, the symmetric bisector plane is the plane that is spanned by e and the average of the facet normals of the two facets adjacent to e. By adding the contributions of all resulting convex bodies using the additivity relationship eq. (9), as outlined in Fig. 6 (c), one obtains the validity of eq. (28) for non-convex triangulated bodies. The sign of the dihedral angle α e ∈ (−π, π] determines if the edge is convex (α e > 0) or concave (α e < 0).

C. Integral Gaussian curvature W 3 (Euler index χ)
As the point-wise Gaussian curvature G 3 on cylinders and flat facets vanishes, only vertices of the triangulation (and their corresponding spherical patches on the parallel body) contribute to W 3 . For both a convex and a non-convex polytope P the point-wise Gaussian curvature G 3 and the integrated Gaussian curvature W 3 can be calculated by the well-known simple sum over angle deficits at surface vertices in eq. (30), derived below, and also given in [71], [72]. The non-convex case is treated by exploiting additivity.
The Gaussian curvature contribution of the vertices v ∈ F 0 is derived by the Gauss-Bonnet-formula where F 2 (v) is the subset of triangles adjacent to vertex v and S denotes the spherical patch on the parallel surface ∂P ǫ . For all ǫ > 0 each spherical cap S ⊂ ∂P ǫ can be uniquely assigned to one vertex v; ∂S its oriented boundary curve and k g the geodesic curvature along ∂S. At the corners of ∂S, the discontinuity of the tangent vectors is characterized by jump angles φ v T (see Fig. 8) which for all ǫ > 0 coincide with the interior angles of the triangle T at v [73], see Fig. 8. The Fig. 9. The Gaussian curvature G 3 at a saddle vertex is obtained by a virtual decomposition of P at v into three polytopes with vertices vµ, µ = 1, 2, 3 using the additivity property of MT. geodesic curvature k g vanishes almost everywhere along ∂S, because ∂S are great circle arcs on the spherical patch and the adjacent cylindrical patch and are thus geodesics; hence the integral ∂S k g ds vanishes.
As a consequence, S G 3 dA is constant for all ǫ > 0. Equation (29) therefore yields a definition and an explicit formula for W 3 (P ) as a sum of the local contribution w 3 (P, v) at a vertex v (30) At a concave vertex v, a polytope P can always be decomposed into three separate bodies (one of vanishing volume) that have convex vertices in lieu of v. It is easy to validate eq. (30) at concave vertices by using the additivity relation in eq. (9) (see Fig. 9). corresponds to the center of mass of P multiplied with its volume (assuming P is homogeneously filled with material of constant density.) The components of this vector may be computed by transforming the volume integral into a surface integral using Gauss' theorem with functions f i that satisfy div f i = x i . The vector-valued auxiliary function f i can be chosen for each index i independently and the index i denotes the index of W 1,0 0 , which is evaluated. For the particular choice of f i given in table III, this can be explicitly written as with k as listed in table III. (k is not a summation index.) |T | is the surface area of T . The I T in eq. (32) are tensorial integrals over the individually parametrized triangles with the three vertices v µ , µ = 1, 2, 3 The components of the auxiliary functions (f i ) k are selected entries of the tensor I T or zero. I T can be written in terms of the triangle vertices v µ and triangle centers C T of T as The remaining integrals W 1,0 ν with ν = 1, 2, 3 are evaluated similarly to the integrals W ν W 2,0 ν (see below). The integrals W 0,1 ν (ν = 1, 2, 3) involving surface normals vanish for arbitrary bodies (with closed bounding surfaces).
The volume integral W 2,0 0 (P ) can be computed in a similar way as W 1,0 0 (P ). Using the components ij of the tensor may be expressed as Again, the index k is not a summation index but rather the index specified in table III. This derivation applies equally to convex and non-convex polytopes P . The computation of W 0,2 1 results in a simple sum of integrals over triangular facets, resulting in the formulae in Tab. II, both for convex and non-convex bodies.
The tensor W 0,2 2 is calculated by a parallel body construction, first demonstrated for convex bodies. Consider a convex polytope P , and the corresponding parallel body P ǫ . The integral over the parallel surface is split up into integrals over flat facets, cylindrical patches and spherical patches. Out of these only the cylindrical edge segments contribute, for the same reasons as for the scalar measure W 2 . The remaining contribution is calculated for ǫ → 0 using the following representation for the normal vectors on the cylindrical patches. Given an edge e with facet normals n T and n T ′ of the adjacent triangles. One obtains (also representing a special case of eq. (48)) To compute the integral on the right hand side we define the orthogonal unit vectorsê = e/|e|,n = (n eT + n e T ′ )/|n eT + n e T ′ | andṅ =ê ×n. For a given edge, n(ϑ) can be written as n = cos ϑn + sin ϑṅ. In this basis, the integral over n 2 evaluates to (1/2) (α e + sin α e )n 2 + (α e − sin α e )ṅ 2 , see Fig. 6. This yields the formula in Tab. II. The validity of this formula for non-convex bodies follows from the same additivity arguments as for W 2 .
G. Curvature-weighted surface integrals W 2,0 2 and W 2,0 3 The mean and Gaussian curvature weighted surface integrals W 2,0 2 and W 2,0 3 over position vectors can be evaluated as the limit ǫ → 0 of the parallel body construction, for convex bodies. The validity for non-convex shapes follows from the analogous construction as for W 2 and W 3 (see tab. II).

H. Open bodies, labeled domains and Minkowski maps
The analysis presented so far has been derived for compact bodies in Ê 3 with a closed bounding surface -and inherits strong robustness from its integral nature. For some analyses, the requirement of closed bodies is too stringent. For example, experimental datasets of percolating or periodic structures, both of which extend infinitely through space, always represent finite subsets of the structure with components that traverse the dataset boundaries. Similarly, an analysis of a periodic model may be restricted to a translational unit cell, see Fig. 10. Furthermore, a local MT analysis, termed a Minkowski map [6], [13], can be useful to quantify variations throughout the sample. For a Minkowski map, a grid is superposed on the body K, and the MT are computed separately for each grid domain L. Such Minkowski maps can be useful to analyze spatial heterogeneity of anisotropy or orientation at the length scale given by the size of L. In these situations, the MT are computed for the subset L∩∂K of the whole bounding surface that is contained in a box L. In general, L∩∂K is not a closed surface even if ∂K is.
It is evidently possible to take the subset K ∩ L of K and consider ∂(K ∩ L) as the bounding surface. However this introduces bounding surface patches (e.g. solid/void interfaces if K is a porous medium) that are not part of the bounding surface ∂K of K. For physical analyses one may want to avoid such boundary effects, i.e. not consider the contributions of these additional bounding surface patches. This motivates the introduction of MF and MT for open bodies, i.e. bodies without a closed bounding surface (see Fig. 10).
In lieu of an attempt to define MF and MT for open bodies, we define a domain-wise analysis of MF and MT. Consider a decomposition of the surface ∂P of a triangulated body P into m domains, or patches, D σ (with σ = 1, . . . , m) such that and consider these domains labeled by labels σ = 1, . . . , m.
Triangles are uniquely assigned to a label, but edges and vertices of the triangulation can be shared between several domains. Specifically, the domains D σ could represent a decomposition of P into patches contained within the grid domains of a three-dimensional lattice (See Fig. 11).
Contributions of the facets can be uniquely assigned to D σ , but the contribution of edges and vertices on the boundaries of a patch D σ needs to be divided between σ and the label of the adjacent domain σ ′ (See Fig. 11).
For W r,s 2 , the contribution of the dihedral angle at an edge is equally divided between the labels of the two adjacent triangles (Note that this is naturally taken account of by the use of oriented edges in the doubly-connected edge list, discussed above). For W r,s 3 the division of the contribution of the interior vertex angles to the integral Gaussian curvature measures W r,s 3 is less straight-forward. An intuitive way, that is also consistent with global integration over all labeled domains, is provided by the use of label factors. The label factor f Dσ (v) of domain D σ at vertex v is defined as where F 2 (D σ ) is the set of all triangles labeled with label σ. Hence f Dσ (v) is the sum of angles at v of those triangles adjacent to v and are labeled σ divided by sum of these angles of all adjacent triangles. It is easy to see, that for a vertex with m adjacent labels and W 3 (P ) = v∈F0 w 3 (v).
For the volume tensor W 2,0 0 a label-wise analysis is only well-defined if the body K is subdivided (and not only the bounding surface ∂K).
(a) (b) Fig. 10. (a) A compact (but not convex) body K that corresponds to a translational unit cell of a periodic body and a triangulation of it. Its bounding surface consists in the translational unit as well as the flat "end caps" that seal it. (b) A surface portion representing a translational unit and its triangulation. The body K is the volume to one side of the surface and forms a connected periodic body. The surface and its triangulation extend beyond the translational unit cell indicating that the surface (and the triangulation) is periodic.

III. ANISOTROPY MEASURES
Based on MT, robust measures of anisotropy can be defined that are sufficiently sensitive to capture subtle anisotropy effects and that are applicable to a wide range of microstructures. The usefulness and versatility of this approach is demonstrated by two examples representing different types of structures -a cellular partition and a network structure.
A rank-2 tensor is defined to be isotropic if and only if it is proportional to the unit tensor Q, i.e. its eigenvalues are all equal. Deviations from isotropy are measured by the anisotropy index β r,s ν , which is the ratio of extremal eigenvalues of the tensor W r,s ν . For example, let ξ µ (|ξ 1 | ≤ |ξ 2 | ≤ |ξ 3 |) be the eigenvalues of W 0,2 1 , then the anisotropy index is By definition, this quantity is dimensionless, continuous and rotation invariant. The value of 1 indicates perfect isotropy, and smaller values indicate anisotropy. For anisotropic bodies, it is sometimes also useful to consider γ 0,2 1 = |ξ 2 /ξ 3 |. These quantities can be easily interpreted for the translation invariant tensors W 0,2 1 and W 0,2 2 . We can write W 0,2 1 equivalently to eq. (22) as the second moment of the distribution of normal vectors (with the density ρ 1 (n)) on the unit sphere Ë 2 as W 0,2 where the ρ 1 (n ′ ) = ∂K δ(n − n ′ ) dA. That is, ρ 1 is areaweighted density of normal vectors. It is easy to see that an uniform distribution on Ë 2 is equivalent to an isotropic tensor. For example, if K is a sphere, then ρ 1 (n) is constant and β 0,2 1 = 1 as expected. For the rectangular box [0, a x ]×[0, a y ]× [0, a z ], the function ρ 1 (n) is concentrated at delta peaks, The resulting anisotropy measure is a z /a x for a x ≥ a y ≥ a z . It is instructive to express the second translation-invariant MT W 0,2 2 by a distribution of normals and curvatures. The density gives the sum of the area of all surface patches that have normal direction n ′ and mean curvature G ′ 2 .
If the function ρ 2 can be written as a product ρ 2 (n, G 2 ) =ρ 2 (G 2 ) ρ 1 (n), the anisotropy characteristics β 0,2 1 and β 0,2 2 defined as the ratio of the smallest to the largest eigenvalue of W 0,2 2 are identical. In this sense, β 0,2 2 provides a higher order anisotropy measure that quantifies anisotropy of the curvature distribution.

A. Alignment of Actin Biopolymer networks under shear
Biopolymer networks made of actin or collagen fibers are important structural elements in biological tissue that act as a scaffolds and provide stiffness and mechanical stability [75]- [77]. Of current interest is the relationship between fiber arrangement and alignment on the one hand side and elastic or visco-elastic properties on the other. This relationship can be probed by shear-experiments with confocal microscopy providing real-space structural data [74]. We now demonstrate that the degree of alignment and of structural anisotropy of the fiber network is well-captured by a MT analysis.
The data sets analyzed here represent actin fiber networks reconstituted from rabbit actin biopolymer networks with actin concentration of 1.2mg/ml cross-linked with filamin A. These are imaged using confocal microscopy for different shear deformations, see the explanation in Fig. 12. The data sets are the same as those analyzed in [74]. The gray-scale data set is converted into a binary data set with 1 corresponding to actin and 0 corresponding to the surrounding fluid by standard threshold segmentation with threshold I c 6 . The medial axis of the 1 phase is computed using distance-ordered homotopic 6 The threshold Ic is chosen such that only the brightest and hence thickest fibers are retained. For a given segmentation threshold Ic the integrated intensity of all voxels of the fluid phase is Λ = ( I(p)) −1 * I(p) where I(p) is the intensity of voxel p in the original intensity data set, the sum over all voxels of the data set and * the sum over all voxels of the fluid phase, i.e. those voxels that are set to 0 by the segmentation process. The values of Ic chosen here correspond to 0.95 < Λ < 0.99. thinning [78], [79] and is used as the one-voxel thick line representation of the actin fiber network. Conversion to a triangulated representation is obtained by using the Marching Cubes algorithm [80]. For more details of the analysis of biopolymers see ref. [81].
Typically only a subset, or observation window, of the structure is available for analysis. Therefore we assume that the network is homogeneous and a sufficiently large but finite subset is accessible which is assumed to represent the entire sample. The derived measures β 0,2 1 quantify the intrinsic anisotropy, i.e. their values do not depend on the size, aspect ratio or position of the observation window (for sufficiently large windows). Figure 12 shows 1/β 0,2 1 and 1/γ 0,2 1 evaluated on the whole network (that consists essentially in a single component; only the outer boundary layers of the confocal data are clipped). It shows that the distribution of normal directions of the fiber bounding surface becomes less isotropic with increasing shear, indicating alignment of the fibers. The angle between the eigenvector to the minimal eigenvalue ξ 1 (corresponding approximately to the dominant tangent direction) and the direction of shear decreases to 0, indicating the alignment of the fibers with the direction of shear. This is commensurate with the results published in [74], that extracted a distribution of tangent directions and used these to quantify alignment.
The eigenvalue ratios of the translation-covariant tensors W 2,0 ν (and of the tensor of inertia) capture different aspects of the anisotropy of a shape compared to the translation invariant tensors W 0,2 ν , see also section I-A. However, usefulness of the translation-covariant tensors depends on whether or not a natural definition of the origin is available for the system. For example, for the analysis of cellular shapes one may choose the center of mass W 1,0 0 /W 0 or the corresponding curvature centroid W 1,0 ν /W ν as the origin. Especially for percolating or periodic bodies, for which the analysis is always restricted to a finite window of observation, the choice of origin is often not naturally determined. An additional problem for such structures is that the measures β 2,0 ν derived from translation covariant tensors W 2,0 ν , as opposed to the translation-invariant measures β 0,2 ν , crucially depend on the shape and size of the window of observation.
The analysis of alignment of biopolymer networks illustrates the potential of the MT approach for structure characterization of cellular and porous materials, and demonstrates its applicability to voxelized experimental data. The MT approach can shed light on systems with a similar spatial structure that exhibit subtle anisotropy effects, such as fibrous biological materials [58], porous materials [4], and metal foams [82].

B. Anisotropy of free-volume cells of random bead Packs
Granular media are a physical system where geometry plays an essential role in determining the physical properties, such as flow or packing properties. The geometric structure to be characterized is substantially different from the above example and consists in an ensemble of grains. These grains do not fill space without gaps, and each grain is associated with its surrounding region of space, often through a construction of the Voronoi diagram.
The distributions of volumes of the Voronoi cells of disordered jammed sphere packs with packing fractions from 0.55 (random loose packing) to around 0.64 (random close packing, RCP) have attracted interest for the study of granular systems [83], motivated by a possible statistical mechanics description of granular systems [84], [85]. Aste et al. [83] used the volume distribution of Voronoi cells to estimate configurational entropy in static packings.
Experimentally, RCP is the maximal packing fraction observed in disordered monodisperse sphere packings. However, for non-spherical particles this limiting packing fraction has a higher value: anisotropic bodies such as M&M chocolate candies pack more tightly in the disordered phase than spherical particles, reaching packing fractions of φ ≈ 0.71 for spheroids and φ ≈ 0.735 for general ellipsoids [86], [87]. Higher order shape measures of the Voronoi cells, such as their anisotropy, have only recently been investigated, demonstrating characteristic anisotropy distributions in Voronoi diagrams of monodisperse spheres [1]. It has been conjectured that ellipsoids are, due to their anisotropic shape, making more efficient use of the anisotropic accessible volume in a static packing and hence pack denser than spheres [1].
Experimental and simulated random close-packed configurations of approximately 10 5 −10 6 monodisperse spheres were generated using three different protocols [88]- [90]. The center coordinates of all beads (in the experimental data sets) were extracted by X-ray microtomography and a deconvolution technique [91], [92]. The Voronoi cells of the sphere centers are computed using the program qhull [93]. For a point p of a set P of points, its Voronoi cell is the convex polytope that contains all points of Ê 3 closer to p than to any other point in P [94]. Figure 13 shows a subset of a experimental bead sample with packing fraction 0.58 on the left panel. The wireframe illustrates the Voronoi diagram of these spheres. On the right hand side, spheres were replaced by ellipsoids with the same anisotropy measure β 2,0 0 and the half-axes a 1 = a 2 ≤ a 3 . a 3 is aligned along the eigendirection of the largest eigenvalue of W 2,0 0 of the corresponding Voronoi cell. The use of MT, and their eigenvalues and eigendirections, hence provides an efficient way of "fitting" ellipsoids to given convex cells.
This substitution of spheres by ellipsoidal particles (allowing slight overlap) leads to packings with φ ≈ 0.72, consistent with experimental and numerical observations for these kinds of packings [87].
This example has demonstrated the ability of MT to quantitatively characterize the subtle, but important, anisotropy effects in packings of granular particles. A quantitative analysis of this system, and of the relevance of anisotropy for an interpretation of maximally obtainable packing fractions in random granular media, is presented in ref. [1].
Voronoi volumes in granular matter can be connected to the configurational entropy [95]; the MT analysis goes beyond the description of volume alone. MT can be used to test local or global isotropy assumptions [1]. This is important for systems with inherent anisotropy, such as dense granular flows under gravity [96] or shear-zones in granular matter [97]. Furthermore, MT (of rank 4) have been proven useful for the precise determination of a crystallization onset in granular systems [7]. This demonstrates the relevance of advanced shape characterization for the study of granular matter.
The characterization of cellular shapes is also relevant in various other contexts, e.g. for the relationship between structure and dynamics in glass-forming liquids [98] or for packing entropy of the hard micellar cores in supramolecular micellar materials [99], and for the understanding of physical properties of foams, both dynamic [100] and static [41], [101].

IV. CONCLUSION
In this article, we have described linear-time algorithms to compute Minkowski tensors of three-dimensional triangulated bodies. We have shown that robust measures of intrinsic anisotropy can be derived from these tensors. These have been used to quantify the anisotropy of an experimental model system for granular matter and of confocal microscopy images of sheared biopolymer networks.
The availability of this algorithm is a prerequisite for applying MT analysis more widely for the quantitative description of anisotropic structural features in physical systems. Further applications of this algorithms in studies concerning anisotropy are discussed in refs. [2], [4] ranging from the analysis of Turing patterns, foams and cellular systems to defect detection in crystals. Rank-4 tensors have already been used to detect crystallization in jammed packings [7].
More and more 3D real space become available in many domains of science. This paper provides a efficient method to extract basic morphological information from experimental 3D data. For the design of optimal spatial structures (see for example refs. [102], [103]) the understanding of the interplay between morphology and physical properties is essential.

SOFTWARE
Software to compute the Minkowski tensors is freely available at our website, www.theorie1.physik.fau.de/karambola ACKNOWLEDGMENT We are grateful to T. Aste, G. Delanay, M. Saadatfar, T. Senden and M. Schröter for the bead pack data, to J. Liu and D.A. Weitz for the actin biopolymer data. We thank C. Marlow (Sigatoka, Fiji Islands) for permission to reproduce his painting "White Spirits", A. Boyde (Imperial College) for the trabecular bone image and M. Saadatfar (Australian National University) for the metal foam image in Fig. 1. We thank Julia Hörrmann and Michael Klatt for their critical comments on the manuscript. GEST, DH and KM acknowledge support by the German research foundation (DFG) through the research group 'Geometry and Physics of Spatial Random Systems' under grants SCHR1148/3-1, ME1361/12-1 and HU1874/2-1.

A. Specific examples
For some simple shapes the rank-2 MT can be calculated analytically by using explicit surface parametrizations and expressions for surface normals and principal curvatures. Specifically for a sphere of radius R centered at the origin, one obtains and, for ν = 1, 2, 3 and r + s = 2, with the unit tensor Q.
For a convex polytope P , we write F µ (P ) for the set of µ-dimensional faces of P , µ = 0, 1, 2, that is, F 0 (P ) is the set of vertices, F 1 (P ) is the set of (non-oriented) edges, and F 2 (P ) is the set of faces. If F ∈ F µ (P ), then we denote by n(P, F ) the set of exterior unit normal vectors of P at F , which is a (2 − µ)-dimensional subset of the unit sphere Ë 2 . Then we obtain, as a special case of general formulas in section I-B, W r,s ν (P ) = 1 3 F ∈F3−ν (P ) F x r H 3−ν ( dx) n(P,F ) n s H ν−1 ( dn) (48) with ν = 1, 2, 3. For the notation see the main text at eq. 17.
For a cuboidal box of dimensions a x × a y × a z aligned with the coordinate axes and centered at the origin eq. (48) yields W 0 = a x a y a z , W 1 = 2 3 (a x a y + a y a z + a z a x ), W 2 = π 3 (a x + a y + a z ), W 3 = 4π 3 , W 1,0 ν = 0 and all MT of rank two are diagonal matrices with the following entries where {i, j, k} = {x, y, z} and permutations thereof.
For an ellipsoid given by (x/l x ) 2 + (y/l y ) 2 + (z/l z ) 2 = 1 the surface integrals all result in elliptic integrals and cannot be expressed in closed form. However, the scalar MF W 0 is W 0 = 4π 3 l x l y l z . The MT W 2,0 0 is diagonal with where {i, j, k} is {x, y, z} and permutations thereof. The integration of all other tensors is easily numerically obtained by use of the ellipsoid parametrization x(u, v) = {l x cos(u) sin(v), l y sin(u) sin(v), l z cos(v)} which yields explicit expressions for the metric tensor of the ellipsoidal surface, the normal vector, and the mean and Gaussian curvatures. These are readily integrated numerically. Fig. 14  as function of lz/lx. In contrast to the above four tensors, this ratio depends strongly on the value of the intermediate radius ly. In particular, for lz = 0 the eigenvalue ratio only becomes zero if the intermediate radius is also ly = 0. For the maximal ly = 1 the eigenvalue ratio converges to 0.5 for lz/lx → 0. minimal to maximal eigenvalue ratio of the MT of rank two of ellipsoids with l x = 1 and 1 ≥ l y ≥ l z .