Local and global energies for shape analysis in medical imaging

In a previous contribution, a new Riemannian shape space, named TPS space, was introduced to perform statistics on shape data. This space was endowed with a Riemannian metric and a flat connection, with torsion, compatible with the given metric. This connection allows the definition of a Parallel Transport of the deformation compatible with the three‐fold decomposition in spherical, deviatoric, and non‐affine components. Such a parallel transport also conserves the Γ‐energy, strictly related to the total elastic strain energy stored by the body in the original deformation. A new approach is here presented in order to calculate the bending energy on the body alone (body bending energy) and to restrict it exclusively within physical boundaries of objects involved in the deformation analysis. The novelty of this new procedure resides in the fact that we propose a new metric to be preserved during the TPS direct transport. This allows transporting the shape change more coherently with the mechanical meaning of the deformation. The geometry of the TPS space is then further discussed in order to better represent the relationship between the Γ‐energy, the strain energy, and the so‐called bending energy densities.

the use of diffeomorphisms, which are more suitable to describe soft tissues changes. [7][8][9] On the other hand, these descriptions are very sophisticated but do not incorporate in their formulation some synthetic properties of the decomposition of deformation in some significant aspects. In this paper, we focus on shapes represented by landmarks that are considered anatomically homologous in different configurations. Very often, in medical applications, the goal of shape analysis is not the assessment of differences between shapes; rather, the main aim is to find differences between deformations. This distinction is crucial and inevitably leads to the definition and quantification of deformation occurring between two shapes, eg, a source and a target. 4,10,11 In this paper, we focus on the left ventricle (LV) cardiac revolution represented by 2D shapes deriving from projections on the coronal plane of the whole LV 3D reconstruction coming from 3D speckle tracking echocardiography (3DSTE). In particular, we present a theoretical framework that suggests an efficient method for calculating the body bending energy (BE) and the Γ-energy observed in a generic deformation. In the present paper, the formulation of the thin plate spline (TPS) space 12,13 is further developed in order to better represent the relationship between the Γ-energy, the strain energy (SE), and the so-called bending-energy densities. The bending energy is then calculated on the body only (body-bending energy) in order to restrict it exclusively within the physical boundaries of objects involved in the deformation analysis. The utility of these quantities is illustrated by an application of the TPS direct transport 12,13 to two LV belonging to one healthy and one pathological individual affected by hypertrophic cardiomyopathy. The use of only two individuals here is purely didactic thus without clinical significance and exclusively aimed at showing the important implications of energy density conservation instead of absolute energy. For computational reasons, we present here a formulation applicable only to 2D cases. However, a 3D generalization can be achieved by further exploring the potential of integration techniques on the surface boundaries.

DEFORMATIONS IN SIZE AND SHAPE SPACE
In landmark-based geometric morphometrics, configurations are represented by k × m matrices, where k is the number of landmarks and m is the dimension of the Euclidean ambient space  m , identified with R m and endowed with an orthonormal basis e i , … , e m fixed once for all. The positions x ∈  m of k points, called landmarks, define a configuration of the body, which can be represented as a k × m matrix X = (x 1 , … , x k ) T , where x i = (x i1 , … , x im ) T are the coordinates of the ith landmark with respect to e i , … , e m . 14 When configurations are centred, a strict representation is given by (k − 1) × m matrices. The space of each centred configuration, the so-called centred configuration space  k m can be identified with the linear space of the (k − 1) × m matrices. On each point X ∈  k m , representing a configuration, the tangent space  X  k m is defined. A vector V ∈  X  k m represents a small deformation of the configuration X and can be represented, in turn, by a (k − 1) × m matrix. We define the centred configuration X C as the k × m matrix where C = I k − 1 k 1 k 1 T k , I k is the k × k identity matrix and 1 k is a k × 1 column of ones. We define the Helmertized configuration X H as the (k − 1) × m matrix where H is the so called Helmert sub-matrix. The jth row of the Helmert sub-matrix H is given by and so the jth row consists of h j repeated j times, followed by jh j and then k − j − 1 zeros, j = 1, … , k − 1. One can switch from one parametrization to the other by using the properties and then If one wants to filter out every rigid motion then rotations are removed, thus obtaining the size-and-shape space SΣ k m =  k m ∕SO(m) as a quotient space with respect to the rotation group. If one wants to filter out the scaling too, then each form is scaled to unit size, thus obtaining the shape space Σ k m . This is the reason upon which the geometry of the shape space depends on the choice of a size measure. This choice is not obvious.
In order to describe and compare different deformations, it is very important to decompose each deformation by projecting it on different subspaces of  X  k m , each one characterizing specific features of the deformation.
In order to describe the deformation, it is imperative to introduce some important subsets of  X ( k m ) introduced in Varano et al. 13 The first is the set  u of the uniform deformations of X It describes the overall deformation of the reference configuration.
There are other two subspaces that are, in turn, subspaces of  u : the spherical (homothetic) deformations of X  sph ⊂  u and the small rotations of X  rot ⊂  u While each one of the above mentioned subspaces has a direct characterization, there are two subspaces, namely, the set  nu of the non-uniform deformations of X, and the set  nsph of the non-spherical deformations of X, which admit only indirect definitions as orthogonal complements. In particular In such a way, the following decompositions hold where ⊕ g is the orthogonal sum with respect to the metric tensor g. In fact, the orthogonality is not a concept given a priori, but it is defined by the metric tensor. We want to stress that g depends on the considered point X then the orthogonality depends on X.
By defining the deviatoric component as the uniform non-spherical component we obtain, in turn, the following orthogonal decomposition, w.r.t. g, on each configuration X (2.8) It is easy to understand that the above decomposition can be interpreted as a splitting of the (small) deformation in a component that changes only orientation, size, shape.
Then, given a small deformation V X of the configuration X, the rotational component V rot must be filtered out by using a proper procedure of alignment (either ordinary procrustes analysis [OPA] or modified ordinary procrustes analysis [MOPA] see below).
If one performs the quotient with respect to the rotation group, for obtaining the size-and-shape space SΣ k m =  k m ∕SO(m), the above mentioned decomposition can be interpreted as a decomposition in vertical subspace ( rot ) and horizontal subspace (the orthogonal complement  sph ⊕ g  dev ⊕ g  nu ). The horizontal subspace on the configuration X represents (is isometric to) the tangent space to the size-and-shape space SΣ k m at the form [X] S , where [X] S is the equivalence class given by the action of SO(m) on X.
It is important to stress that the above decomposition, in general, relies on the choice of g. In the next sections, we will introduce, as g, the so called TPS metric. 12

Thin plate spline
By following the notation of Varano et al 12 and Dryden and Mardia, 14 we summarize the construction of the TPS and we refer to previous studies 12,14,15 for further details. In the Euclidean space  m , the m-tuple of interpolating TPS is a function represented by the pair (A X , W X ), where: A X is a linear transformation of  m , represented by a (m × m) matrix and W X is a (k × m) matrix. Given a point x ∈  m , and a centred configuration X ∈  k m , we have where s(x) = ( (x − x 1 ), … , (x − x k )) T a is (k × 1) matrix, x i ∈ X is the position of the ith landmark, and Given an undeformed configuration X, and a deformed configuration X ′ , we can apply Equation (3.10) landmark-wise, yielding to There are 2k interpolation constraints in Equation (3.11), and we introduce m × (m + 1) more constraints on W in order to uncouple the affine and non affine parts: For a given pair (X, Y), there exists a unique set of parameters for the pair (A X , W X ) that solve the problem (3.11), constrained with (3.12): are a m × (k − 1) and a (k − 1) × (k − 1) matrices, respectively, which only depend on the source configuration X. Finally, the target X ′ can be represented as the deformation of X 14) The centred coordinates can be recovered simply by pre-multiplying with H T . The matrix (k − 1) × (k − 1) matrix Γ 11X is called the bending energy matrix and represents the nonlinear part of the deformation. As it vanishes on linear deformations, its eigenvectors associated to nonvanishing eigenvalues are only (k−1−m). They are called principal warps and represent the principal modes of deformation. 12,14

Strain energy, bending energy, body bending energy
One of the main features of the TPS method is that it minimizes the so-called BE J, involving the second derivative of the displacements. 15 In particular, given a body, represented in the undeformed configuration by the regular region Ω of R m , we label x the points in Ω and x ′ the points in the deformed configuration Ω t at the instant t. If the configurations are sampled in k landmarks, then Ω can be represented by the k × m matrix X and Ω t by the deformed configuration X. The displacement field is represented by the difference vectors In the landmark case, the displacements experienced by the k landmarks can be collected in the The bending energy is defined as follows: Note that the bending energy is defined as an integral on the whole R m , which does not diverge because the curvature of the TPS, outside Ω, decays quickly to zero. 16 The bending energy can be considered as a pseudo-metric on   k m because it measures the difference between two configurations but it vanishes on the affine deformations.
A nonsingular distance, used in the deformable templates method, 17 is the Dirichlet energy D .
A more physical distance, used in continuum mechanics, which vanishes on the rotational part of the local deformation, is the complete where and are the Lamé elastic moduli, related to the material properties of the body, and is the strain tensor (in the case of small displacements). Here, as we are describing only geometrical features, we assume a simplified expression for the SE In Varano et al, 13 in order to give a physical interpretation of the bending energy, the body-bending energy (BEB) has been introduced, performing the integration only on the body: In general, the BE will be slightly greater than the BEB and the decay = J∕J Ω can be used to quantify this difference. In Varano et al, 13 it has been shown that for bilinear deformations in 2D, in which a rectangle bends in a trapezoid where = 1 e 1 + 1 e 2 are the bending with respect to the two axes, the strain energy is directly proportional to the body-bending energy and the proportionality coefficient depends only on geometrical properties of Ω.
In particular, the following relationship holds: where  p and  are the polar inertia and the area of Ω, respectively. On the other hand, for a uniform deformation, such that x ′ = Ax, and A is a (m × m) matrix, it is clear that the strain energy is directly proportional to the square of the Frobenius norm of (A − I) where I is the (m × m) identity matrix. A general bilinear deformation can be parameterized as follows:

BE matrix and BEB matrix
An interesting and useful feature of the TPS is that the value of the BE J associated to a deformation U of a configuration X is directly proportional to the quadratic form Tr(U T Γ 11X U) For this reason, the matrix B ∶= Γ 11X is also called the BE matrix of X. In fact, one of the most important reasons of the success of the TPS is the possibility to evaluate the BE directly by means of a closed form expression, avoiding the need to discretize the configuration in a huge number of finite elements (triangles in 2D and tetrahedra in 3D). Furthermore, the BE matrix depends only on the undeformed configuration and can be used as pseudo-metric on   k m . We show here that also for the BEB calculation, it is possible to define a symmetric matrix B Ω , holding Then the decay (U) can be calculated by means of the Rayleigh quotient:  In order to build the BEB matrix, we follow, using a different notation, the procedure described in Kent and Mardia. 16 In the following, for the sake of simplicity, we limit ourselves to consider 2D configurations defined by landmarks on the boundary of the body and in the interior (Figure 1).
The whole procedure is based on the following property of the biharmonic function : where o is the Dirac function defined as the generalized function such that As explained in the previous sections, once performed a TPS on two configurations X and X ′ = X + U, we can describe the deformation as: In order to calculate the BEB, we take the second gradient thus obtaining: The BEB can be then calculated by the following integral: where the BEB matrix is then defined as follows: In order to calculate the component of the matrix, we integrate by parts: where n is the outer unit normal vector field to the boundary. Then, by using Equations (3.29) and (3.30), we obtain for the bulk term of the integral where the values of the angles i , taking into account the integration of the Dirac delta on Ω instead of R 2 , depend on the position of the ith landmark with respect to Ω, in particular Figure 1). For the boundary term, we obtain where the p th side of Ω is represented by x p + p , n p is the unit normal to the pth side, * ( p ) is the vector p rotated clockwise by ∕2.We note that q is the number of landmarks that does not lie on the boundary Ω, while (k − p) is the number of landmarks lying on Ω. Finally, we can put together each term It is possible to prove (see Kent and Mardia 16 for the proof) that fixing the landmarks and changing the integration domain, in such a way that Ω → R 2 , then B Ω → 16 B. This can be qualitatively understood by considering that in case the boundary term C Ω vanishes, and for the bulk term C Ω , there are no landmarks belonging to Ω and then i = 2 for each i, and by using the notable property Γ 11X S X Γ 11X = Γ 11X . 12,13 It is important to stress that for Ω ⊂ R 2 , the bulk term C Ω has the closed form representation (3.35), while the boundary term C Ω is calculated solving, numerically, the single variable integrals appearing in Equation (3.37), which involve the derivatives of the biharmonic function. The BEB matrix computation is here presented only for 2D data because the procedure of integration by parts in that case leads to a line integration that can be easily performed. In the 3D case, Equations (3.35) to (3.37) should become much more complicated. However, a 3D generalization can be achieved by exploiting numerical integration techniques on a surface's boundaries.

TPS metric and Γ-energy
In this section, we introduce a modified expression for the TPS metric formulated in Varano et al 12 and discussed in another study. 13 This modification is motivated by two reasons: the introduction of the BEB matrix and the consideration that the best quantity to conserve in a parallel transport of deformation is not the total strain energy but the mean strain energy density. Lets consider two bodies having the same shape but a different size (volume) undergoing the same deformation, for example, a global stretching of 30%, storing different quantities of SE, but the same average SE.
Then, we propose a general representation for the metric tensor: It is important to stress that Γ 21X , B Ω , 2 depend on the source configuration X. Let given two configurations X and X ′ = X + U, by means of the introduced metric, we can define the distance between X and X ′ called Γ-Energy, as follows: We could consider the Γ-energy Γ(U) as the average strain energy (Ũ) evaluated on a more simple deformationŨ characterized by the same uniform componentŨ u = U u of U and a bilinear deformationŨ nu storing the same body bending energy of U nu , ie, such that J Ω (Ũ nu ) = J Ω (U nu ).

Filtering rotations: OPA and MOPA
Once introduced the TPS metric tensor G, a general technique for aligning two forms [X] S and [Y] S is based on the minimization of the distance defined as The aligned configurationŶ is obtained by means of an optimal rotationQ minimizing d.
). According to this definition,Q turns out to be the rotational component of the polar decomposition of Y T GX. When = 0, we obtain the rotational component of Y T X (the classical OPA 14 ). When = 1 then Y T GX = A X . In the latter case, we define the alignment modified OPA (MOPA). 12

TPS direct transport
In Varano et al, 12 a connection named TPS connection was introduced. It has the following properties: 1. It is compatible with the TPS metric. 2. It is compatible with the decomposition given in (2.9). 3. It preserves the Γ-energy. 4. It is independent of the path.
In this section, we reformulate this parallel transport procedure by implementing it with a new feature that takes in consideration the physical boundaries of the bodies involved. We also assume all the configurations to be centred, and represented by the Helmertized landmarks, then, if not otherwise specified, each matrix is a (k − 1) × m matrix. Furthermore, deformation vectors will have a subscript denoting the starting point, that is, the source configuration (for details see Varano 12 ).
Let X and Y be two source configurations and V X and V Y the two associated deformation vectors, given by (3.42) We say that V Y is the parallel transport of a given V X , that is, and the non uniform part W Y of V Y solves the linear systems:  where the (k − 1) × (k − 1 − m) body principal warps matrix E X , collects all the body principal warps of X, while Q X is a suitable (k − 1 − m) × (k − 1 − m) orthogonal matrix (ie, Q T X Q X = I), defined on each configuration X, representing a rotation, or a reflection, of the principal warps. Once chosen a configuration P as a Pole for the space, Q X is estimated minimizing the Euclidean distance ||E X Q  X − E P || between the rotated principal warps of X, and the corresponding basis on the pole P.
The first equation of (3.43) constrains W b to be orthogonal to the affine part, while the second define the isometry in the subspace of the non affine deformations. This last requirement implies the conservation of the total bending energy. The system (3.43) can be written as Which can be re-written as and so The Equation (3.44) characterizes V Y as the parallel transport of V X . The only difference between the original TPS direct transport 12 and that introduced here is in the definition of the (k − 1) × (k − 1 − m) matrix E X , collecting, here, all the body principal warps of X in columns. Then we call it the body principal warps matrix. We can build it as follows: • Perform a TPS analysis on X and find the BEB matrix B Ω (Equation 3.38).
• Perform an eigenvalue analysis on B Ω and obtain B Ω = ΓΛΓ T where Γ is the (k − 1) × (k − 1) matrix containing, in column, the eigenvectors i , and Λ is the diagonal (k − 1) × (k − 1) matrix of the eigenvalues 1 , … , k−1 , ordered by increasing magnitude. The first m eigenvalues will be equal to 0. • Drop the first m columns from Γ, by obtaining the (k − 1) × (k − 1 − m) matrixΓ, containing the body principal warp eigenvectors by column. • Drop the first m rows and the first m columns from Λ, by obtaining the  The same steps must be used to build the body principal warps matrix E Y on the target configuration. All the above mentioned formulations have been implemented in R by using ad hoc built and easy to use R functions available from the authors upon request.

CARDIAC EXAMPLES
We tested our procedure on real data coming from 3D-STE of human LV. In order to illustrate the approaches we describe here, we contrasted two complete cycles relative to one healthy subject and to a patient affected by hypertrophic cardiomyopathy. The use of only two individuals here is purely didactic thus without clinical statistical meaning and exclusively aimed at showing the profound implications of the body energy density conservation instead of total bending energy. These data come from the same research project partially published in previous studies, [18][19][20][21][22][23][24] and have been collected by means of 3D-STE (PST25SX Artida, Toshiba Medical Systems Corp., Tokyo, Japan). The whole LV geometry is reconstructed by starting from a set of six homologous landmarks, manually detected by the operator for each subject under study. The manual detection of the initial set of landmarks is crucial to fulfil a homology principle, so that comparable anatomical structures of different subjects are detected; completely automated approaches may suffer from error of pattern identification depending on the specific algorithm used for reconstruction. The data set generated for each subject by 3D-STE consists of a sequence of 3D configurations, acquired every 45 ms, each constituted by 1297 landmarks, assumed FIGURE 7 TPS deformation grids for the healthy individual computed from the healthy heart cycle' Grand Mean to be homologous, for both the epicardial and endocardial surfaces, positioned along 36 horizontal circles, each comprised of 36 landmarks, plus the apex (Figure 2 left). After the generation of the data set, we define 16 time events according to Piras et al, 18 to be used as homologous times to sample the heart beating (using cubic spline interpolation) for both subjects as shown in Figure 3, and we assemble a trajectory for each subject by extracting from the data set a sequence of 32 configurations corresponding to the homologous time-events. 18 Finally, in order to present our results in 2-dimensions, for each trajectory, we project a coronal slice of the 3D configurations (interpolated at homologous times) on the plane transversal to the LV base, identified by three points on the epicardium: the LV apex, the point lying on the mitral annulus in correspondence on the middle of the interventricular septum and that lying on the opposite lateral side (Figure 2 right); this plane is selected at the first homologous time, ie, at the R peak (=telediastole) and all shapes corresponding to other times were projected on that plane. This allows to better conserve, during projection, the torsional movement of the LV for the duration of the cycle. Our procedure entirely mimics the 2D-speckle tracking (2DSTE) approach except for the lesser time' resolution. However, this does not affect our result's interpretation given the methodological aim of the present contribution.

Bilinear deformations of a vertical cross section of a LV
A simulated dataset (Figure 4 ) was generated by deforming, parametrically, the diastolic state of the healthy individual through bilinear deformations (Equation 3.25). Deformations were built using six parameters, A 11 , A 12 , A 21 , A 22 , 1 , 2 each following six series of thirty randomly chosen values, and successively sorted in ascending order, in a uniform [0, 0.01], 2 ∈ [0, 0.01], respectively. The diastolic state was chosen as the source configuration (=the undeformed shape) while the 30 deformed shapes were the targets. We then computed the BEB matrix on the source and we used it in order to compute the body bending energy and the gamma energy according to (3.33) and (3.41). We als'o computed the strain energy and the bending energy according to (3.20) and (3.16).

Deformations of a set of vertical cross sections of acquired LVs
We used the two cycles in order to proceed with an implemented version of the TPS Direct Transport procedure described in Varano et al. 12 We chose the consensus of the all 32 shapes enlarged five times (in order to generalize the procedure by considering also a large size difference) as the common reference (CR) toward which per-cycle deformations, estimated from per-cycle local references (LOCS), were transported. We chose as LOCS the two consensus of each individual cycle. In order to achieve this, we computed the BEB matrices of CR and LOCS as well as 1 CR, 2 CR, 1 LOCS vector (two values in this case as we have two LOCS), and mu2LOCS vector (two values). We then contrasted the energy densities computed after the TPS direct transport procedure between the common reference and transported data contrasted with the energy densities calculated on original data between LOCS and their proper deformed states. We expect isometry for both BEB density and gamma energy density, while it is not expected for SE density and BE density.

FIGURE 9a
Thin plate spline (TPS) deformation grids for the transported data computed from the chosen common reference toward which deformation were transported FIGURE 9b Continued Figure 5 left shows that the SE is very close to the Γ-energy in the simulated example ( coefficient: 1.14); in particular, at the beginning of the deformations series, their relationship is very close to isometry, while it is not exactly the case for larger deformations. This could be due to the triangulation's discretization during SE calculation. Figure 5 right shows, instead, that the BEB (on y-axis) is much smaller than the BE calculated on the whole R m . Figure 6 (top right) shows that the Γ-Energy density calculated on original data between LOCS and their proper deformed states are perfectly conserved on transported data when computing it from CR. Figure 6 (top left) shows also the non-affine part of the Γ-energy density that can be obtained by multiplying the BEB density by the proper 2 (i.e. that belonging to the corresponding LOC) coefficient for each deformed state. This part also is perfectly recovered by the TPS direct transport. Bottom side of Figure 6 illustrates that neither the BE nor the SE are perfectly recovered. For the former case, in fact, BEB is formally different from the BE. The SE density is not perfectly conserved as, being the real case not produced by parameterized bilinear transformations, it does not equate the Γ-energy density. However, the deviation from isometry is not dramatic especially for the healthy individual (ie, the black circles: Beta coefficient: 1.12). Figures 7 and 8 show the transformation grids of original data relatively to their proper LOC of the healthy and pathological subject respectively. Figure 9 shows the transformation grids of transported data relatively to the chosen CR. Figure 10 shows the utility of a parallel transport. The upper panel shows the PCA space following a GPA: the PC1 explains the majority of variance and is related mainly to the difference between subjects, while the differences between deformations is relegated to the PC2 that explains a very small fraction of total variance. On the other hand, as shown in the FIGURE 10 PCA and deformations associated to PC axes' extremes computed on original data (upper panel) and on transported data (lower panel). Deformation grids shows colour according to log(det(I + ∇u)) coming from the two-dimensional thin plate spline first derivative evaluated within the body. Values less than 0 indicate that, with respect to the source (here the sample consensus in red), the target (here the PC's extremes in black) experiences a reduction in the local area, while values greater than 0 indicate an enlargement bottom panel, after the use of TPS direct transport we can better appreciate the differences in deformation, ie, magnitude, direction and size of the two trajectories. Deformation grids shows colour according to log(det(I + ∇u)) coming from the two-dimensional TPS first derivative evaluated within the body. Negative values indicate that, with respect to the source (here the sample consensus in red), the target (here the PC's extremes in black) experiences a reduction in the local area, while values greater than 0 indicate an enlargement.

DISCUSSION AND CONCLUSIONS
We have shown that restricting the calculation of the BE to the body, instead using it calculated on the whole R m , allows us to tune the original TPS direct transport in order to conserve Γ-energy and BEB densities. The BEB is estimated by building the corresponding BEB matrix that can be thought as a modification of the classic BE matrix but having the advantage that the physical non-affine deformation of the body is now correctly captured within its real boundaries. This implementation of the original formulation of the TPS direct transport makes it particularly suitable in cases where the analysis of deformation is oriented to a mechanical interpretation of the phenomenon under study. We also stress that formulations presented here have been developed for 2D cases but they are extendable to three-dimensions in order to embrace a wider range of examples and to better approach reality for those studies that can exploit modern 3D imaging technologies.

APPENDIX A
In this section, we summarize some basic algebraic operations used in the paper. In particular, let  m and  m be the Euclidean space of dimension m and the related translation space, respectively. We identify  m ≡  m ≡ R m through the choice of an origin o and a orthonormal basis e i , … , e m . Each vector v ∈ R m , can be represented as where i = 1...m and summation convention over repeated indexes is employed. The dot product between two vectors u, v ∈ R m is then defined as A second order tensor is a linear map on  m ; the dot product between two second order tensors A, B ∈ R (m×m) is defined as

A · B = Tr(A T B).
We can define second order tensors by composing two vectors through a tensor product; by definition, the tensor product of a, b ∈ R m is the linear map on  m ; defined by In the same way, the tensor product between three vectors a, b, c ∈ R m is a third order tensor, here defined by means of its action on a fourth vector u ∈ R m :