Minkowski tensors of anisotropic spatial structure

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 triangulation consisting of planar polygonal facets. For a polytope, the tensors W^2,0_ν characterize the distribution of mass if the cell is a solid cell (W^2,0₀), a hollow cell (W^2,0₁), a wire frame (W^2,0₂) and a cell consisting of points at the vertices only (W^2,0₃); in the last two cases, however, this distribution of mass is weighted with certain exterior angles (see figure 3).

Zoom In Zoom Out Reset image size

The tensor of inertia I_P, defined by ${I}_{P_{ij}}=\int _K \left (-\mathbf {x}_i\,\mathbf {x}_j+\delta _{ij}|\mathbf {x}|^2\right ) \,{\mathrm { d}} V$, is a measure of the mass distribution of a (homogeneous) body K, relevant for the relationship between a rotation and the resulting moment. As I_P is not translation-covariant, it is not a linear combination of the MT. However the simple relationship ${I_P}(K)=-W_0^{2,0}(K)+\mathrm {tr}\left (W_0^{2,0}(K)\right )Q$ holds for arbitrary K; as above Q is the unit tensor of rank 2 of $\mathbbm {R}^3$, also called metric tensor. This illustrates that W^2,0₀(K) is a measure of the distribution of mass if K is a homogeneously filled solid, somewhat analogous to the tensor of inertia. Similarly, W^2,0₁(K) characterizes the mass distribution if K is hollow and bounded by an infinitesimally thin surface sheet. This interpretation of W^2,0₀ and W^2,0₁ is valid for bodies K bounded by arbitrary surfaces (not just polytopes).

For a polytope P, the tensor W^2,0₂(P) reduces to a line integral over the edges of the polytope (as the mean curvature vanishes on the flat facets, see also section 2) 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₃(P) reduces to a sum of point contributions, as the Gaussian curvature G₃ of P vanishes except at the vertices of the given polytope P. Again due to the additivity requirement, these vertices need to be weighted with the Gaussian curvature G₃.

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 2.

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 [20, 22, 23, 63].

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 figure 4 for d = 2) of a convex body K in d-dimensional Euclidean space $\mathbbm {R}^d$ to a suitable local construction.

Zoom In Zoom Out Reset image size

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 $\mathbf {x}\in \mathbbm {R}^d$. This defines a continuous map $\mathbf {p}_K(\cdot ):\mathbbm {R}^d\to K$, x↦p_K(x). Then d_K(x): = ∥x − p_K(x)∥ is the distance from x to K and n_K(x): = (x − p_K(x))/d_K(x), for $\mathbf {x}\in \mathbbm {R}^d\setminus K$, is an exterior unit normal of K at the boundary point p_K(x)∈∂K (see figure 4). This construction achieves a definition of nearest points (reminiscent of the Euclidean distance map [64]) and surface normals that is also well defined for points x whose nearest point on K is a sharp corner (where 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 [23] to obtain local support measures. The definition of MF and MT then follows directly as a special case [33, 62]. The intuitive idea underlying the definition of a local parallel set is described in two steps. Firstly, we specify a region $L\subset \mathbbm {R}^d$ and some  > 0. Then we consider all points $\mathbf {x}\in \mathbbm {R}^d\setminus K$ which have distance d_K(x) at most from K and for which p_K(x)∈L. Secondly, 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\subset \mathbbm {S}^{d-1}$ of the unit sphere in $\mathbbm{R}^d$. As an example, consider the polytope P in figure 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 $S\subset \mathbbm {S}^{d-1}$ of normal directions, conveniently combined to $\eta = L\times S \subset \mathbbm {R}^d \times \mathbbm {S}^{d-1}$. This two-step procedure can be merged and slightly extended to the following general definition.

For given  > 0 and $\eta \subset \mathbbm {R}^d\times \mathbbm {S}^{d-1}$, the local parallel set of K defined by

$$\begin{equation} \mathcal{M}_\epsilon (K,\eta):=\{\mathbf{x}\in \mathbbm{R}^d\setminus K:d_K(\mathbf{x})\leqslant \epsilon,(\mathbf{p}_K(\mathbf{x}),\mathbf{n}_K(\mathbf{x}))\in\eta\} \end{equation} \tag{ 13 }$$

contains all points $\mathbf {x}\in \mathbbm {R}^d$ with 0 < d_K(x) ⩽  such that the pair (p_K(x),n_K(x))∈η. The first condition restricts x to the global outer parallel set $(K\uplus B_\epsilon)\setminus K$ (see figure 4(b1)). The second condition restricts $\mathcal {M}_\epsilon (K,\eta )$ to those points $\mathbf {x}\in \mathbbm {R}^d$ of the global outer parallel set, where (p_K(x),n_K(x))∈η.

The volume of this local parallel set of a convex body K is $V_d(\mathcal {M}_\epsilon (K,\eta ))$, where V_d is the d-dimensional Lebesgue-measure. A fundamental result in integral geometry, known as the local Steiner formula [23, 63], states that the map $\epsilon \mapsto V_d(\mathcal {M}_\epsilon (K,\eta ))$ for all  > 0 is a polynomial of degree d, that is

$$\begin{equation} V_d(\mathcal{M}_\epsilon(K,\eta))=\sum_{\nu=0}^{d-1}\epsilon^{d-\nu}\kappa_{d-\nu}\Lambda_\nu(K,\eta), \end{equation} \tag{ 14 }$$

where $\kappa _n:=\pi ^{n/2}/\Gamma (\frac {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 equation (14) to be true for all  > 0, it is crucial that K is convex [65]. Equation (14) is easily confirmed (and evaluated) for a convex polytope P. In this case, the set $(K\uplus B_\epsilon)\backslash P$ can be decomposed in an elementary way into wedges over the faces F of P as indicated by figure 4(b1). Let $\mathcal {F}_\nu (P)$ denote the ν-dimensional faces of the polytope P⁵

$$\begin{equation} \mathbf{n}(P,F):=\{\mathbf{n}_P(\mathbf{x})\in \mathbbm{S}^{d-1} : \mathbf{p}_P(\mathbf{x}) \in \mathrm{relint}\, F , \mathbf{x}\in\mathbbm{R}^d \backslash P \} \end{equation} \tag{ 15 }$$

the set of unit normal vectors assigned to $F\in \mathcal {F}_\nu (P)$ and relint F the relative interior of F (i.e. the interior of F w.r.t. the lowest-dimensional embedding affine hull); see figure 5.

Zoom In Zoom Out Reset image size

The contributions to $V_d(\mathcal {M}_\epsilon (K,\eta ))$ coming from each of these wedges

$$\begin{equation} \mathcal{M}_\epsilon(P,F):=\mathcal{M}_\epsilon(P,(\mathrm{relint}\,F\times \mathbbm{S}^{d-1})) \end{equation} \tag{ 16 }$$

can be calculated by a simple integration known as Fubini's theorem or Cavalieri's principle [66, 67]. This shows that for η = L × S equation (14) holds with

$$\begin{equation} \Lambda_\nu(P,L\times S)=\frac{1}{\omega_{d-\nu}}\sum_{F\in \mathcal{F}_\nu(P)} \int_{F\cap L}\,{\rm d} \mathcal{H}^\nu\int_{\mathbf{n}(P,F)\cap S}\,{\rm d} \mathcal{H}^{d-\nu-1}, \end{equation} \tag{ 17 }$$

where ω_ν = ν κ_ν is the surface measure of the (ν − 1)-dimensional unit sphere and $\,{\mathrm { d}} \mathcal {H}^\nu$ is the Hausdorff measure of dimension ν [68].

Since an arbitrary convex body K can be approximated by polytopes, equation (14) can be derived by a continuity argument. In fact, Λ_ν(K,η) can be expressed as a linear combination of $V_d(\mathcal {M}_{\epsilon _\mu }(K,\eta ))$, for _μ > 0 pairwise different and μ = 1,...,d. One obtains an invertible d × d-matrix equation (with the matrix entries ^d−ν_μκ_d−ν) with ν running from 0,...,d − 1. Therefore, the properties of the local parallel volume $V_d(\mathcal {M}_\epsilon (K,\eta ))$ (in particular, additivity and weak continuity) are also available for Λ_ν(K,η) [63, p 202]. In particular, K↦Λ_ν(K,η) is an additive functional for fixed η. That is, for convex bodies K and K'

$$\begin{equation} \Lambda_\nu(K \cup K',\eta)=\Lambda_\nu(K ,\eta)+\Lambda_\nu(K',\eta)-\Lambda_\nu(K \cap K',\eta). \end{equation} \tag{ 18 }$$

Furthermore, η↦Λ_ν(K,η) is a non-negative measure for fixed K. The latter means that if $\eta _\mu \subset \mathbbm {R}^d\times \mathbbm {S}^{d-1}$, $\mu \in \mathbbm {N}$, is a sequence of mutually disjoint (measurable) sets, then

$$\begin{equation} \Lambda_\nu\left(K,\bigcup_{\mu=1}^\infty\eta_\mu\right)=\sum_{\mu=1}^\infty\Lambda_\nu(K,\eta_\mu). \end{equation} \tag{ 19 }$$

This property is called σ-additivity of Λ_ν(K,·). Weak continuity of Λ_ν means that for every sequence of convex bodies K_μ $(\mu \in \mathbbm {N})$, with K_μ → K and every continuous function $f:\mathbbm {R}^d\times \mathbbm {S}^{d-1}\rightarrow [0,\infty )$ the equation

$$\begin{equation} \lim_{\mu\rightarrow \infty} \int f(\mathbf{x},\mathbf{n}) \, \Lambda_\nu(K_\mu,d(\mathbf{x},\mathbf{n}))\nonumber = \int f(\mathbf{x},\mathbf{n}) \,\Lambda_\nu(K,d(\mathbf{x},\mathbf{n})) \end{equation} \tag{ 20 }$$

holds. Note that this does not imply Λ_ν(K_μ,η) → Λ_ν(K,η) as μ → ∞ for all measurable sets $\eta \subset \mathbbm {R}^d\times \mathbbm {S}^{d-1}$

In particular, Λ_ν(K,·) can be used to integrate functions over $\mathbbm {R}^d\times \mathbbm {S}^{d-1}$. It is plausible that Λ_ν(K,·) is concentrated on the normal bundle $N(K):=\{(\mathbf {p}_K(\mathbf {x}),\mathbf {n}_K(\mathbf {x}))\in \partial K\times \mathbbm {S}^{d-1}:\mathbf {x}\in \mathbbm {R}^d\setminus K\}$. In other words, N(K) consists of all $(\mathbf {x},\mathbf {n})\in \partial K\times \mathbbm {S}^{d-1}$ such that n is an exterior unit normal vector of K at x. The measures Λ_ν(K,·), ν = 0,...,d − 1, are called support measures and are determined as coefficient measures of the Steiner formula, equation (14). They are local versions of the classical MF W_ν, since $\Lambda _\nu (K,\mathbbm {R}^d\times \mathbbm {S}^{d-1})=V_\nu (K)\propto W _{d-\nu }(K)$.

If K is sufficiently smooth, then

$$\begin{equation} \Lambda_{\nu}(K,\eta)=\frac{{{d-1} \choose {\nu}}}{\omega_{d-\nu}} \int_{\partial K} \mathbf{1}\left\lbrace(\mathbf{x},\mathbf{n}_K(\mathbf{x}))\in\eta \right\rbrace G_{d-\nu}(\mathbf{x})\, \,{\rm d} A, \end{equation} \tag{ 21 }$$

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 evaluates to one if its argument is true and 0 otherwise. For general dimensions and ν ⩽ d − 1, equation (21) holds for all sufficiently smooth 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

$$\begin{equation} \Phi^{r,s}_\nu(K):= \frac{ 1}{r!s!}\frac{\omega_{d-\nu}}{\omega_{d-\nu+s}} \int_{\mathbbm{R}^d\times\mathbbm{S}^{d-1}} \mathbf{x}^r\mathbf{n}^s\,\Lambda_\nu\left(K,d(\mathbf{x},\mathbf{n})\right) , \end{equation} \tag{ 22 }$$

hence we obtain Φ^r,s_ν(K) by integrating the tensorial function x^rn^s with respect to the measure Λ_ν(K,·) over $N(K)\subset \mathbbm {R}^d\times \mathbbm {S}^{d-1}$. If K is a polytope, for d = 3 this yields equation (A.3), up to a different normalization. If K is smooth, we obtain equation (5), up to a different normalization and indexing scheme, i.e.

$$\begin{eqnarray} & &\Phi^{r,s}_\nu(K)=\frac{{{d-1} \choose {\nu}}}{r!s!\omega_{d-\nu+s}}\int_{\partial K}\mathbf{x}^r\mathbf{n}^s G_{d-\nu}(\mathbf{x})\, \,{\rm d} A,\nonumber\\ & &\Phi^{r,0}_d(K)=\frac{1}{r!}\int_{K} \mathbf{x}^r \mathrm{d}V. \end{eqnarray} \tag{ 23 }$$

The notation Φ_μ,r,s or Φ^r,s_μ for the MT in equation (22) is preferred in some of the mathematical literature and differs from the notation W^r,s_ν in equations (4) and (5) only by a different indexing scheme and a different normalization. In $\mathbbm {R}^3$, i.e. for d = 3, the functionals Φ^r,s_μ and W^r,s_ν are related by

$$\begin{equation} W_\nu^{r,s}(K)=C\,\Phi_{d-\nu}^{r,s}(K)\quad \mathrm{with}\; C:=\frac{r!s!\omega_{\nu+s}}{d{{d-1} \choose {\nu-1}}}, \end{equation} \tag{ 24 }$$

for ν = 1,...,d, and

$$\begin{equation} W_0^{r,0}(K)=r!\,\Phi_{d}^{r,0}(K). \end{equation} \tag{ 25 }$$

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 equation (22).

Since the theory of support measures is well developed [23, 63], 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

$$\begin{eqnarray} &&\fl\int_{\mathbbm{R}^d\times\mathbbm{S}^{d-1}} \mathbf{x}\mathbf{n}\,\Lambda_\nu\left(K\uplus\mathbf{t},d(\mathbf{x},\mathbf{n })\vphantom{\int}\right)=\int_{\mathbbm{R}^d\times\mathbbm{S}^{d-1}} (\mathbf{x}+\mathbf{t})\mathbf{n}\,\Lambda_\nu\left(K,d(\mathbf{x},\mathbf{n} )\vphantom{\int}\right)\nonumber\\ &&=\int_{\mathbbm{R}^d\times\mathbbm{S}^{d-1}} \mathbf{x}\mathbf{n}\, \Lambda_\nu\left(K,d(\mathbf{x},\mathbf{n})\vphantom{\int}\right)+\mathbf{t}\int_{\mathbbm{R}^d\times\mathbbm{S}^{d-1}} \mathbf{n}\,\Lambda_\nu\left(K,d(\mathbf{x},\mathbf{n})\vphantom{\int}\right). \end{eqnarray} \tag{ 26 }$$

It is a basic property of the measures Λ_ν(K,·) that they are centred at the origin in the sense that $\int \mathbf {n}\,\Lambda _\nu (K,d(\mathbf {x},\mathbf {n}))=0$ [34], 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 equation (22) to measurable subsets of $\mathbbm {R}^d\times \mathbbm {S}^{d-1}$.

The calculation of the volume of a polytope P can be transformed into a surface integral by the Gauss theorem, see equation (11) and [70]. With div x = div(x,y,z)^t = 3 one obtains

$$\begin{equation}\fl W_0(K)=\int_P \,{\rm d} V = \frac{1}{3}\int_P \mathrm{div}\, \mathbf{x} \,{\rm d} V =\frac{1}{3}\int_{\partial P} \langle \mathbf{x}, \,{\rm d}\mathbf{A}\rangle =\mathrm{tr}\, \frac{1}{3}\int_{\partial P} \mathbf{x}\mathbf{n}\,\mathrm{d}A= \mathrm{tr}\, W_{1}^{1,1},\end{equation} \tag{ 27 }$$

where  dA = n dA denotes the oriented infinitesimal area element and 〈·,·〉 the scalar product.

The surface integral is a sum over triangles and is easily evaluated yielding the formulae in table 2. This result is independent whether P is convex or not. The surface area W₁(P) of ∂P is simply the sum of triangle areas.

Expressing W₂(P) as the limit of vanishing parallel distance of W₂(P) of the parallel body P, W₂(P) = lim_→0W₂(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 $\int \int _S \frac {1}{\epsilon } \, \epsilon ^2\sin \theta \,\,{\mathrm { d}}\varphi \,\,{\mathrm { d}}\theta$ 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

$$\begin{equation} W_2(P) = \frac{1}{3}\lim_{\epsilon\rightarrow 0} \sum_{{\bf e}\in\mathcal{F}_1} \frac{1}{2}\int_0^{|{\bf e}|} \,{\rm d} l \int_{-\alpha_{\bf e}/2}^{\alpha_{\bf e}/2} \frac{1}{2\epsilon} \epsilon\,{\rm d}\vartheta =\sum_{{\bf e}\in\mathcal{F}_1} |{\bf e}|\,\frac{\alpha_{\bf e}}{12} ,\end{equation} \tag{ 28 }$$

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 figure 6. Note that $\mathcal {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.

Equation (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 (i.e. −π < α_e < 0), see figure 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 equation (9), as outlined in figure 6(c), one obtains the validity of equation (28) for non-convex triangulated bodies. The sign of the dihedral angle α_e∈(− π,π] determines whether the edge is convex (α_e > 0) or concave (α_e < 0).

Zoom In Zoom Out Reset image size

As the point-wise Gaussian curvature G₃ on cylinders and flat facets vanishes, only vertices of the triangulation (and their corresponding spherical patches on the parallel body) contribute to W₃. For both a convex and a non-convex polytope P the point-wise Gaussian curvature G₃ and the integrated Gaussian curvature W₃ can be calculated by the well-known simple sum over angle deficits at surface vertices in equation (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 $\mathbf {v}\in \mathcal {F}_0$ is derived by the Gauss–Bonnet formula

$$\begin{equation} \int_S G_3 \,{\rm d} A = 2\pi-\sum_{T\in\mathcal{F}_2(\mathbf{v})} \phi_T^{\mathbf{v}}-\int_{\partial S} k_{\mathrm{g}} \,{\rm d} s, \end{equation} \tag{ 29 }$$

where $\mathcal {F}_2(\mathbf {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 figure 8) which for all  > 0 coincide with the interior angles of the triangle T at v [73], see figure 8. The 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 $\int _{\partial S}k_{\mathrm {g}} \,{\mathrm { d}} s$ vanishes.

Zoom In Zoom Out Reset image size

As a consequence, $\int _{S} G_3 \,{\mathrm { d}} A$ is constant for all  > 0. Equation (29) therefore yields a definition and an explicit formula for W₃(P) as a sum of the local contribution w₃(P,v) at a vertex v

$$\begin{equation} W_3(P)=\frac{1}{3}\sum_{\mathbf{v}\in\mathcal{F}_0} w_3(P,\mathbf{v})=\frac{1}{3}\sum_{\mathbf{v}\in\mathcal{F}_0} \left( 2\pi-\overset{}{\underset{T\in\mathcal{F}_2(\mathbf{v})}{\sum}}\phi_{T}^{\mathbf{v}}\right). \end{equation} \tag{ 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 equation (30) at concave vertices by using the additivity relation in equation (9), see figure 9.

Zoom In Zoom Out Reset image size

The Minkowski vector W^1,0₀ corresponds to the centre 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

$$\begin{equation} \left(W_0^{1,0}(P)\right)_{i} = \int_{P}{} \mathbf{x}_i \,{\rm d} V= \int_{\partial P}{} \langle \mathbf{f}_i , \,{\rm d} \mathbf{A} \rangle \end{equation} \tag{ 31 }$$

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₀, which is evaluated. For the particular choice of f_i given in table 3, this can be explicitly written as

$$\begin{equation} (W_0^{1,0}(K))_{i} = \overset{}{\underset{T\in\mathcal{F}_2}{\sum}} \int_{T}{} \mathbf{x}_i \mathbf{x}_k (\mathbf{n}_{T})_{k} \,{\rm d} A = \overset{}{\underset{T\in\mathcal{F}_2}{\sum}} (I_{T})_{ik} (\mathbf{n}_{T})_{k} \end{equation} \tag{ 32 }$$

with k as listed in table 3. (k is not a summation index.) |T| is the surface area of T. The I_T in equation (32) are tensorial integrals over the individually parametrized triangles with the three vertices v_μ, μ = 1,2,3

$$\begin{equation} I_{T} = 2\vert T\vert \int_0^1 \,{\rm d} a \int_0^{1-a} \,{\rm d} b \left[\mathbf{v}_{1}+a(\mathbf{v}_{2}-\mathbf{v}_{1})+b(\mathbf{v}_{3}-\mathbf{v }_{1})\right]^2. \end{equation} \tag{ 33 }$$

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 centres C_T of T as

$$\begin{equation} I_{T} = 2\vert T\vert \left(\frac{9}{24} \mathbf{C}_{T}^2 + \frac{1}{24} \overset{3}{\underset{\mu=1}{\sum}} \mathbf{v}_{\mu}^2 \right). \end{equation} \tag{ 34 }$$

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₀(P) can be computed in a similar way as W^1,0₀(P). Using

$$\begin{equation} J_{T} = 2\vert T\vert \int_0^1 \,{\rm d} a \int_0^{1-a}\,\,{\rm d} b \left[\mathbf{v}_{1}+a(\mathbf{v}_{2}-\mathbf{v}_{1})+b(\mathbf{v}_{3}-\mathbf{v }_{1})\right]^3, \end{equation} \tag{ 35 }$$

the components ij of the tensor may be expressed as

$$\begin{equation} \left(W_0^{2,0} (K) \right)_{ij} = \overset{}{\underset{T\in\mathcal{F}_2}{\sum}} (J_{T})_{ijk} (\mathbf{n}_{T})_{k}. \end{equation} \tag{ 36 }$$

Again, the index k is not a summation index but rather the index specified in table 3. This derivation applies equally to convex and non-convex polytopes P.

The computation of W^0,2₁ results in a simple sum of integrals over triangular facets, resulting in the formulae in table 2, both for convex and non-convex bodies.

The tensor W^0,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₂. 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 equation (A.3),

$$\begin{equation} W_2^{0,2}(K) = \frac{1}{12}\sum_{\mathbf{e}\in \mathcal{F}_1} \vert\mathbf{e}\vert\int_{-\alpha_{\mathbf{e}}/2}^{\alpha_{\mathbf{e}}/2}\mathbf{n}^2 \,{\rm d}\vartheta. \end{equation} \tag{ 37 }$$

To compute the integral on the right-hand side we define the orthogonal unit vectors $\hat {\mathbf {e}}=\mathbf {e}/\vert \mathbf {e}\vert$, $\ddot {\mathbf {n}} =(\mathbf {n}_{\mathbf {e}_T}+\mathbf {n}_{\mathbf {e}_{T'}})/\vert \mathbf {n}_{\mathbf {e}_T}+\mathbf {n}_{\mathbf {e}_{T'}}\vert$ and $\dot {\mathbf {n}}=\hat {\mathbf {e}}\times \ddot {\mathbf {n}}$. For a given edge, n(υ) can be written as $\mathbf {n}=\cos \vartheta \,\ddot {\mathbf {n}}\kern-1pt+\kern-1pt\sin \vartheta \,\dot {\mathbf {n}}$. In this basis, the integral over n² evaluates to $\frac12\left ((\alpha _{\mathbf {e}}\kern-1pt+\kern-1pt\sin \alpha _{\mathbf {e}})\ddot {\mathbf {n}}^2\kern-1pt+(\alpha _{\mathbf {e}}\kern-1pt-\kern-1pt\sin \alpha _{\mathbf {e}})\dot {\mathbf {n}}^2\right)$, see figure 6. This yields the formula in table 2. The validity of this formula for non-convex bodies follows from the same additivity arguments as for W₂.

The mean and Gaussian curvature weighted surface integrals W^2,0₂ and W^2,0₃ 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 for W₂ and W₃ (see table 2).

The analysis presented so far has been derived for compact bodies in $\mathbbm {R}^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 data sets of percolating or periodic structures, both of which extend infinitely through space, always represent finite subsets of the structure with components that traverse the data set boundaries. Similarly, an analysis of a periodic model may be restricted to a translational unit cell, see figure 10. Furthermore, a local MT analysis, termed a Minkowski map [6, 15], 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 analyse 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.

Zoom In Zoom Out Reset image size

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 figure 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 $\sigma =1,\ldots ,m$) such that

$$\begin{equation} \partial P= \bigcup_{\sigma=1}^{m}D_\sigma, \end{equation} \tag{ 38 }$$

and consider these domains labelled by labels $\sigma =1,\ldots , 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 3D lattice (see figure 11).

Zoom In Zoom Out Reset image size

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 figure 11).

For W^r,s₂, 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 into account by the use of oriented edges in the DCEL, discussed above). For W^r,s₃ the division of the contribution of the interior vertex angles to the integral Gaussian curvature measures W^r,s₃ is less straightforward. An intuitive way, that is also consistent with global integration over all labelled domains, is provided by the use of label factors. The label factor f_Dσ(v) of domain D_σ at vertex v is defined as

$$\begin{equation} f_{D_\sigma}(\mathbf{v}):=\frac{ \sum\nolimits_{T\in \mathcal{F}_2(\mathbf{v})\cap\mathcal{F}_2 (D_\sigma)}{}\phi_{T}^{\mathbf{v}} }{ \sum_{T\in \mathcal{F}_2(\mathbf{v}) }{}\phi_{T}^{\mathbf{v}}}, \end{equation} \tag{ 39 }$$

where $\mathcal {F}_2(D_\sigma )$ is the set of all triangles labelled with label σ. Hence f_Dσ(v) is the sum of angles at v of those triangles adjacent to v and are labelled σ divided by the sum of these angles of all adjacent triangles. It is easy to see that

$$\begin{equation} w_3(\mathbf{v})=\sum_{\sigma=0}^{m} f_{D_\sigma}(\mathbf{v}) w_3(\mathbf{v}) \end{equation} \tag{ 40 }$$

for a vertex with m adjacent labels and $W_3(P)=\sum _{\mathbf {v}\in \mathcal {F}_0}w_3(\mathbf {v})$.

For the volume tensor W^2,0₀ a label-wise analysis is only well defined if the body K is subdivided (and not only the bounding surface ∂K).

A fully functional implementation of the algorithms represented in sections 2.1–2.8 is provided as online supplementary data (available from stacks.iop.org/NJP/15/083028/mmedia) to this paper, and also made available through the internet at www.theorie1.physik.fau.de/karambola, under a GNU General Public License.

The implementation is a straightforward realization of the formulae in the rightmost column of table 2 into ANSI-C code. A simple data structure is used to store a triangulation of a surface, as a set of points and a list of facets (specifically triangles); the data structure allows us to iterate over all vertices, edges or facets by simple loops, for example over all edges 'for e in $\mathcal {F}_1$', and it allows the extraction of neighbours, for example of all (one or two) triangles that are adjacent to a given edge e. With that data structure in place, the sums in table 2 simply become for loops, that are all linear in the number of edges, facets or vertices.

Biopolymer networks made of actin or collagen fibres 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 fibre 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 fibre network is well captured by a MT analysis.

The data sets analysed here represent actin fibre networks reconstituted from rabbit actin biopolymer networks with actin concentration of 1.2 mg ml⁻¹ cross-linked with filamin A. These are imaged using confocal microscopy for different shear deformations, see the explanation in figure 12. The data sets are the same as those analysed in [74]. The grey-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.⁸ The medial axis of the 1 phase is computed using distance-ordered homotopic thinning [78, 79] and is used as the one-voxel thick line representation of the actin fibre network. Conversion to a triangulated representation is obtained by using the Marching Cubes algorithm [80]. For more details of the analysis of biopolymers see [81].

Zoom In Zoom Out Reset image size

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₁ 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₁ and 1/γ^0,2₁ 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 fibre bounding surface becomes less isotropic with increasing shear, indicating alignment of the fibres. The angle between the eigenvector to the minimal eigenvalue ξ₁ (corresponding approximately to the dominant tangent direction) and the direction of shear decreases to 0, indicating the alignment of the fibres 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 1.1. 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 centre of mass W^1,0₀/W₀ 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 exhibits subtle anisotropy effects, such as fibrous biological materials [82], porous materials [4] and metal foams [83].

Granular media represent a system where geometry plays an essential role in determining its physical properties, such as flow or packing properties. The geometric structure to be characterized is substantially different from the above example. It consists in an assembly of (disjoint) grains that have, at most, mutual point contacts.

A commonly used way to associate each grain with its corresponding region of space is the construction of a 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) have attracted interest for the study of granular systems [84], motivated by a possible statistical mechanics description of granular systems [85, 86]. For instance, Aste et al [84] used the volume distribution of Voronoi cells to estimate configurational entropy in static packings, and Zhao et al [87] to quantify spatial correlations in disc packings.

Here we illustrate how MT can be used to characterize the shape, rather than simply the volume, of the grains' Voronoi cells, in the spirit of [1].

Figure 13 shows a subset of an experimental data set of static disordered monodisperse jammed spheres with packing fraction 0.58 on the left panel (for details see [1]) The wire-frame illustrates the Voronoi diagram of these spheres⁹. On the right hand side, spheres were replaced by ellipsoids that have the same eigenvalue ratio of the Minkowski tensor W^2,0₀ as their Voronoi cells, W^2,0₀(ellipsoid) = W^2,0₀(Voronoi cell), implying in particular the same value of the anisotropy measure β^2,0₀. Their half-axes are aligned with the eigendirections of the tensor W^2,0₀, evaluated for the corresponding Voronoi cell. (We have here chosen the origin to coincide with the center of mass W^1,0₀/W₀ of the Voronoi cell.) The use of MT, and their eigenvalues and eigendirections, hence provides an efficient way of 'fitting' ellipsoids to given convex cells and hence a means to quantify their anisotropy or elongation. This in turn provides sensitive tools to quantify shape and structure of both amorphous and ordered particulate systems and packings.

Zoom In Zoom Out Reset image size

Quantitative analyses of the Voronoi cell shapes of disordered sphere configurations, in various phases, have been given in [1, 2, 7, 90], illustrating the breadth of the Minkowski tensor approach and in particular its usefulness to capture the onset of crystallization. An application to ordered sphere packings has been given in preliminary form in the conference proceeding [91].

The scope of MT for the analysis of granular material is, however, not restricted to the detection of local crystalline domains. Rather, as the following analysis illustrates, these shape measures also allow for a quantitative description of the local structure in amorphous assemblies, discriminating sharply between different types of amorphous geometries.

Figure 14 shows a study where the simple anisotropy measure β^2,0₀ derived from the MT reveals a distinct difference between the different phases that hard sphere systems can adopt: the plot shows the typical shape, quantified by β^2,0₀, of a Voronoi cell of a given volume, expressed as the local packing fraction ϕ_l: = W₀(grain)/W₀(Voronoi cell). The plot contains data for equilibrium hard sphere systems in the fluid and in the ordered phase (the same data as used in [2]) and of jammed static sphere packs (the same data sets mentioned above and used in figure 6 of [1]). For each sphere configuration, the Voronoi cells of all particles are computed and their local packing fractions ϕ_l. Then all cells are classified by their value of ϕ_l, and the average cell shape 〈β^2,0₀〉(ϕ_l) is determined by averaging β^2,0₀ over all cells of a data set that have, up to a discretization interval Δϕ_l, the same value of ϕ_l. The curves 〈β^2,0₀〉(ϕ_l) describe Voronoi cell anisotropy as a function of the local packing fraction ϕ_l.

Zoom In Zoom Out Reset image size

Figure 14 elucidates how MT can help discern some of the morphological differences between the different hard sphere phases. Firstly, this analysis of the average cell anisotropy clearly discerns the distinct geometries between the static jammed packings and the (also disordered) equilibrium fluid configurations and, more expectantly, the differences from the equilibrium ordered phase. Computing 〈β^2,0₀(ϕ_l)〉 provides a signature of the origin of the data sets, clearly discerning the structure of equilibrium hard spheres and of the static jammed packings.

Secondly, the data for 〈β^2,0₀〉(ϕ_l) for all six jammed static configurations collapse to a single curve that is approximately linear. These data sets comprise different global packing fractions and preparation protocols, including tomography data of dry acrylic beads [95] and glass beads settled against a fluid current [95], as well as Lubachevsky–Stilinger simulations of frictionless particles without gravity [96] and discrete element method simulations of spheres with friction and gravity [97]. This suggests a universality of the jammed static disordered spheres: the typical shape (quantified by β^2,0₀) of a Voronoi cell of a given local packing fraction ϕ_l is the same for all packings, regardless of protocol and of the global packing fraction ϕ_g. Note that this is in stark contrast to the case of the equilibrium fluid or ordered phase, where particularly small cells in a globally denser packing are more isotropic with larger β^2,0₀ than the same size cell in a looser packing. This result suggests that the global packing fraction of a jammed static disordered sphere packing is the result of combining typical 'building blocks' with a given local ϕ_l in different proportions. This observation may go some way towards clarifying why random close packing of spherical particles appears to be largely protocol-independent.

Figure 14 has provided a test case where MT concisely discern the differences between two types of amorphous structures, in a quantitative fashion with an intuitive geometric real-space interpretation. The characterization of amorphous 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, including rheological [100] and static [101], properties, and the evolution under ageing or coarsening [43, 102]. MT can be applied to these systems in an analogous fashion to the analysis of this section.