Topological Approach for Material Structure Analyses in Terms of R 2 Orientation Distribution Function

: The application of solid mechanics theory for material behavior faces the discrete nature of modern or biological material. Despite the developed methods of homogenization, there are deviations between simulated and experiments results. The reason is homogenization, which mathe-matically involves a type of interpolation. The situation gets worse for complex structured materials. On the other hand, a topological approach can help in such analysis, but such an approach has computational costs. At the same time, increasing modern computational capabilities remove this barrier. This study is focused on building a method to analyze material structure in a topological sense. The orientation distribution function was used to describe the structure of the material. The plane case was investigated. Quadratic and biquadratic forms of interpolant were investigated. The persistent homology approach was used for topology analysis. For this purpose, a persistence diagram for quadratic and biquadratic forms was found and analyzed. In this study, it is shown how scaling the origin point cloud inﬂuences H 1 points in the persistence diagram. It was assumed that the topology of the biquadratic form can be understood as a superposition of quadratic forms. Quantitative estimates are given for ellipticity and H 1 points. A dataset of micro photos was processed using the proposed method. Furthermore, the supply criteria for the interpolation choice in quadratic or biquadratic forms was formulated.


Introduction
Solid mechanics assume continuum space. Of course, most materials differ from this state. For these purposes, different types of materials' symmetry are considered (e.g., tetragonal, orthorhombic, etc.). So, quantifying the material structure has become an important task, because the concrete type of symmetry should be found for the current material. Moreover, even if the type of symmetry is known, the corresponding mechanical properties should be found. This task has become crucial for composite and biological materials. So, one of the popular approaches is mathematical modelling [1][2][3]. In these cases, the microstructure's mechanics are simulated directly, and then some general relation is investigated. Then, for real models, this relation is usually used [4]. In other words, the homogenization approach is widespread. The weak points of such methods are computational labor and the idealization of the microstructure. Of course, in models, some idealized geometry is used, and real objects deviate from the ideal, which causes difficulty. For example, bone structure has a complicated architecture, and in different points of the bone different structures appear. The most common measurements have a histological origin, and, in fact, they have mean or relative characteristic [5,6]. So, some clusterization by tissue type of the object is performed. Then, isotropic [7,8] or orthotropic properties [9] and continuum mechanics are used. It is obvious that such an approach suffers from large spread. However, the received values of mechanical properties vary a lot, and modern approaches are focused on personalization of the tissue models, such as computed tomography [10] and digital volume correlation [11]. Computed tomography allows not only the simulation of the object behavior but also analyses of the structure. So, the alternative is the orientation distribution function [12][13][14]. Such an approach compares the structure with some distribution and then the distribution can be interpolated. Then, the interpolation is used for measuring the structure. Furthermore, the weak point of the approach is the type of interpolation. Thus, quadratic form interpolation is well known, but it does not suit all distributions. The background of the quadratic interpolation is a three (two in R 2 ) orthogonal materials symmetry planes hypothesis. Or, in other words, there is a hypothesis of orthotropic materials. The objectivity of the hypothesis is in question. Contrariwise, the topological approach for data analysis has become widespread [15][16][17]. Furthermore, such an approach allows for quantifying the shape of the structure and, consequently, for understanding the material symmetry. Despite this, there is a lack of research dedicated to material structure analyses via the topological approach. So, this study is focused on topological analyses of the material structure. The R 2 case was investigated, and the orientation's distribution function was used for the material structure description. Quadratic and biquadratic interpolations were considered.

Orientation Distribution Function
Let us consider the macroscopic scalar property of the material ϕ ( Figure 1). Usually, the ϕ depends on the material distribution in the sample. The distribution depends on position x and orientation vector n. So, ϕ can be defined by mapping from unit sphere to reals: The scalar-valued function f should satisfy the condition: The invariance of f leads to existing function h: So, h can be expanded in the following Fourier series [18][19][20]: where where F is the first two tensor spherical harmonics and they form the orthogonal basis of the functions.
The three coefficients in the basis can be determined from f (n) via the integrals [21]: In mechanics, it is widespread that the second order tensor is used to describe the microstructural anisotropy in orthotropic materials. There is also the so-called mean intercept length (MIL) tensor [22] or fabric tensor [23,24] ( Figure 2). The MIL tensor is defined as: Furthermore, the fabric tensor is defined as: Such an approximation can be useful to determine the elastic tensor: where φ-porosity, H-fabric tensor, and ρ-scale factor depending on the materials.

Application of Orientation Distribution Function
Obviously, the interpolant of the orientation distribution function can be applied for fabric distribution, e.g., for micro photos. Previously, the authors developed a method for the analysis of the quality of the derma's recovery. For this purpose, collagen distribution in the area was analyzed from the orthotropic media point of view [25,26].
The method's idea is to compare the micro photo and the vector field. The vector field consists of fabric tensor eigenvectors normalized by the degree of anisotropy, where the degree of anisotropy is an eigenvalues ratio: In this case, the quality of the derma's recovery can be estimated by the following equation: In other words, the relative content of the fabric tensor with a low degree of anisotropy is measured. Such an approach assumes that the equality of all eigenvalues leads to isotropic distribution. This was shown by Cowin [27], who developed a general representation of an elastic tensor as a function of the solid volume fraction and of the invariants of the fabric tensor. Furthermore, it was shown that eigenvectors of the fabric tensor are orthotropic axes.
So, the question regards the quality of MIL interpolation by the fabric tensor. Furthermore, as the distribution is badly interpolated by the constant, there are distributions that are badly interpolated by the fabric tensor. On the other hand, four order interpolation is harder to analyze. This is why it was decided to develop a method for distribution analyses.

Topological Approach
Let us define a function ζ : T → R (12) where T-topological space. In this case, T a can be defined as: With the obvious property: Such inclusions produce a map of homologous groups Considering the sequence: We receive: This leads to homeomorphism: This homeomorphism can endure topological classes from T ai to T aj . Topological classes can appear (birth), disappear (death), or unite. The image of homeomorphism keeps this information; denote it: The Betti number can be calculated as rank: The simplicial complex is used, so let us consider the simplicial complex K and denote ζ as: We receive: and In these cases, the birth of class [c] of homology H p in K i can be defined as The death of class [c] of homology H p in K i can be defined as The birth and death of class [c] are illustrated in Figure 3. So, let us describe some technical notes. We consider a point cloud: In the first step, a distance matrix should be built: where d is a Euclidean distance between points: The received distance matrix can be understood as an adjacency matrix for a constant value of d ij . So, each unique value of d ij generates a finite graph. Let us collect unique values in ascending order in a vector d. In this case, the adjacency matrix for d K -balls can be built as: where θ-Heaviside step function. A point cloud and graphs for unique distances are shown in Figure 4. An adjacency matrix for the shown point cloud is: where the vector of unique distances is: The incidence matrix can be built using the adjacency matrix and the unique distances vector. Each unique distance produces a set of lines ordered by the minimal point index. An incidence matrix for the example is: Let us consider the birth and death of class [c] for the point cloud in the example from [28]. Initially, there are 5 classes (each point in Figure 4); d 1 produces the death of two classes (list of classes: point 1, point 5, and the complex of points 2, 3, and 4 in Figure 4). Only one class is left for d 2 (the complex of all connected points in Figure 4). Then, there appears for d 3 new classes, and the list of classes is: the complex of connected points 1, 2, 4, and 5 and the complex of connected points 2, 3, and 4. In the next step (diameter d 4 ), the complex of connected points 1, 2, 4, and 5 breaks into two complexes: the connected points 1, 2, and 4 and the complex of connected points 1, 4, and 5. Each complex of connected points is a 1-dimensional complex. Starting with a d 5 diameter, 1-dimensional complexes disappear. In the 6 th step, complex 2, 3, and 4 became closed, which means that the simplex is covered with d 5 spheres. So, the 1-dimensional complex is dead. Complexes 1, 2, and 4 and 1, 4, and 5 are dead on d 6 and d 7 , respectively.
The incidence matrix can be built for 1-dimensional complexes:  5). The 1-dimensional complex's death is described by the corresponding diameter value. So, knowing the incidence matrices, B0 and B1 information about the birth and death of 0-dimensional and 1-dimensional complexes can be obtained.
To find a birth of a 1-dimensional complex, a cycle should be found in the B0 graph. To find a death of a 1-dimensional complex, a diameter overlapping the whole cycle should be found. This problem can be decomposed by overlapping the simplex. For this purpose, the following equation was used: where x i is a vector of the simplex vertices, D is the diameter of the overlapping, and x is a point of the circles' intersection.
The common way to represent a complex's birth and death is through a persistence barcode and a persistence diagram. So, the persistence barcode for 0-dimension complexes of the example are shown in Figure 5a. Each interval starts from the birth of the complex and ends in a point corresponding to the death diameter [29,30]. In Figure 5b, the persistence barcode for 1-dimension complexes of the example is shown. Understanding the diameters of the birth and death as points in R 2 leads us to a persistence diagram (Figure 5c). It is obvious that points have a sense if they are higher than the bisector; otherwise, they immediately were born and died (if a point is on the bisector) or they died before birth (if a point is under the bisector) [31,32]. Scaling the point cloud leads to scaling the 0-dimension persistence barcode and shifting the 1-dimension persistence barcode. In Figure 5c, the scaled points are shown by circles with black contours. So, the general algorithm of the construction of B0 and B1 matrixes is presented in pseudocode below Algorithm 1.

Algorithm 1: B0 and B1 matrixes construction
Input: x, vector of dimension N Output: B0 and B1 matrixes for each xi for j in range [i, N] d(i,j) = distance(x(i), x(j)) by Equation (28) end for end for Where B0 and B1 are arrays to store the incidence matrixes for H 0 and H 1 , respectively; d is the distance matrix; dUnice is the vector with unique distance values; V, E, and F are temporary variables to store the number of vertices, edges, and faces for a certain distance; and for the cycle store graphs cycles, D stores the overlapping distances for the cycles.
In the case of numerical calculations, some noisy results appear. In Figure 5c, a bisector is highlighted. If a point appears on a bisector, it means that the birth and death of the complex was immediate. If a point appears under a bisector, it means that the complex died before birth. Such points should be filtered. To filter the persistence diagram, bottleneck distance was used:

Research Design
Consider the 2nd order interpolation (quadratic form) of MIL in the following form: The biquadratic form for the 4th order interpolation [33] was used in the following equation: Assuming M 2 in diagonal form (space is oriented by the eigenvectors of M 2 ), two parameters of M 2 and three parameters of M 3 were varied. According to (2) and (6), only antisymmetric and closed curves were inspected [34,35].
The point clouds of curves were obtained by varying the parameters in (36) and (37). Schematically, the received curves are shown in Figure 6. All curves were scaled in a square [−0.5, 0.5] × [−0.5, 0.5]. A persistence diagram was built for each curve.
Persistence diagrams were filtered near the bisector using bottleneck distance (35) with relative value 1%.
where s-scale factor, I-identity matrix, and y-new coordinates. The Euclidean distance (28) will be changed as follows: So, all unique distances (31) will be multiplied by s; this means that all H 0 points will be shifted in a persistence diagram in the death direction.
The birth of graph cycles depends only on unique distances. So, all birth points will be multiplied by scale s. Consider the overlapping Equation (34) in scaled axes: This leads to multiplying by scale s the overlapping diameter. Because the multipliers of the birth and death dimeters are equal, the aspect ratio will be constant.

Lemma 2.
Birth and death values when scaled by λ (0 < λ < 1) in one direction point cloud in R 2 cannot be higher than the origin values.
Proof of Lemma 2. Let us consider changes in distances in scaled dimension. Define the origin distance as d 0 : Then, estimate the distance in scaled dimension d λ : By considering the d λ as a function of λ, monotony can be shown for 0 < λ ≤ 1: The decline of the distance values in the scaled space means that the distance values in the adjacency matrix (29) cannot be higher than in the origin space. That means that the death values of H 0 and the birth values of H 1 cannot increase. Moreover, overlapping distance values (34) cannot be higher than in origin space, too. Furthermore, that means that death values of H 1 cannot increase.

Numerical Example
Then, 10 µm tissue sections were cut with a microtome (Thermo Scientific, Waltham, MA, USA, HM325) to obtain 5 µm thick cuts. These were dehydrated and stained either according to the Mallory protocol or with hematoxylin/eosin.
The histological samples were documented on a Carl Zeiss Axio Imager 2 microscope with magnification of 100×. The image size was about 4164 × 3120 with 150 dpi. The images were taken from a dataset [26]. On the micro photos, the skin biopsy was presented. The biopsy data include the epidermis, the dermis, and the subcutaneous panniculus carnosus muscle. The Cartesian grid was built for images. The cell size was 400 × 400 pixels with 150 dpi. The binarization threshold was calculated by Otsu's method for a blue channel.
To estimate the direction of the collagen orientation, the histological images were aligned with a 2-D Cartesian grid, and the MIL was calculated for each element [36]. Previously, the analysis was obtained through the aspect ratio of fitting ellipse eigenvalues. In this study, the persistence diagram was built for received MILs. MIL distribution was calculated by the following equation: where L-is an intercept length in n direction, ω is a binary function, 1 is for material, and 0 is for no material.

Results and Discussion
Numerical investigation of a quadratic form's persistent diagram shows linear dependence between the death-birth aspect ratio and ellipticity (see the blue line in Figure 7). Let us define the slope value as a death-birth aspect ratio. Linear regression was used for interpolation (see the orange line in Figure 7): The received parameters of linear regression are shown in Table 1.  According to Lemma 1, the slope values are constant for scaling the whole data. This fact allows for the use of the linear regression results for quantifying the ellipticity by points in a persistent diagram.
Numerical investigation of a biquadratic form's persistent diagram shows a significant difference in the case of cross-like curves (e.g., the green curve in Figure 6). In this case, the second H 1 point appears in the persistent diagram. Meanwhile, all quadratic forms follow by one H 1 point.
Some typical point clouds can be highlighted (see Figure 8c) for the received MIL distribution (47) of derma micro photos. Original (see Figure 8a) and binarized (see Figure 8b) cells from micro photos are shown, too. Quadratic interpolation was performed (36), and in Figure 8c, corresponding eigenvectors are shown. It should be noted that received point clouds are not all well-interpolated by quadratic form (e.g., the cyan point cloud in Figure 8). This points to the fact that for analysis, eigenvalues and eigenvectors are used, and some received interpolants are crucially mistaken. However, on the other hand, the persistence diagram for the cyan point cloud in Figure 8 shows two H 1 points (see cyan points in Figure 9). The same is true for cross-like biquadratic curves.  On the other hand, persistence diagrams for biquadratic forms allowed for identifying the areas of curve shape persistence. Numerical results show that flattering shifts the point on the persistent diagram to the right. Moreover, a second H 1 point appears when a new direction of flattering appears (see Figure 9). The area under the bisector is highlighted by a grey color. For example, MIL distribution (the cyan point cloud in Figure 8c) can be implied as a superposition of two ellipses with different aspect ratios of semi axes. In a persistent diagram, these are expressed by two cyan points in Figure 9a.
In contrast, the green point cloud from Figure 8c is the closest to the circle, and the corresponding point on the persistent diagram (the green point in Figure 9) is the closest to the ideal (the left side of the red area). Next, the MILs (the red and blue point clouds in Figure 8c) are flattered in one direction, and these lead to shifting to the right of points in the persistent diagram (the red and blue points in Figure 9a). The values of the slope for lines separating the area by flattering degree are presented in Table 1. Slope values for the H 1 point of the cyan point cloud were 2.5266 and 1.4194, respectively. Using inverse regression, the calculated ellipticity was 0.3609 and 0.1490, respectively. The interpolations by the two ellipses with the following flattering are shown in Figure 9b.
So, the persistent diagram of the concrete MIL point cloud can be understood as decomposition by the number of ellipses. In this case, the number of H 1 points describes the number of composite ellipses, and slope values allow for quantifying the flattering of composite ellipses. The distance from origin quantifies the scale of the composite ellipses. On the one hand, according to this information, MIL distribution can be classified; on the other hand, the required interpolant can be used. Additionally, the number of H 1 points can provide information about the materials' symmetry. So, one H 1 point (or elliptic interpolation) indicates orthotropic material; in contrast, more than one H 1 point raises questions about the orthotropic hypothesis.
Focusing on micro photos analysis from a topological point of view, some items should be highlighted. The number of H 1 points in the persistent diagram indicates the suitable interpolation order. This indicates the symmetry of the MIL distribution; consequently, this also indicates the material symmetry and the orthotropic aspect ratio. With regard to the derma analysis, this information allows for the estimation of the quality of healing (by Equation (11)). Cells with two H 1 points can be analyzed by elliptic decomposition or in terms of biquadratic interpolant. Numerical results showed that two H 1 points appear in cells where collagen distribution deviates from normal. Of course, these facts need to be additionally investigated in proper ways.

Conclusions
In this article, a persistent diagram was built for a quadratic form set and for an antisymmetric biquadratic form set. It was found that the number of H 1 points determines the shape of the curve. Meanwhile, each H 1 point can be understood as an ellipse. The ratio of coordinates is a slope coefficient, which determines the flattering ratio of the ellipse. Linear regression was found for ellipticity and slope values. It was shown that the points shift with the constant slope value describes the scale factor of the origin point cloud. More than one H 1 point determines the biquadratic origin of the data. Furthermore, slope values can be used for interpolation by the number of ellipses. The proposed method was applied on a micro photo dataset. The material was analyzed in terms of the mean intercept length distribution. Distinctive point clouds were presented and discussed in terms of the topology structure. Slope values for applying the proposed method were presented.
Of course, most presented results are numerical, but they are a background for more certain analytical research. More detailed analyses of the topological structure of the antisymmetric biquadratic form can improve our understanding of corresponding persistent diagrams and our understanding of the symmetry of materials (in terms of distribution function), as a consequence. So, the aforementioned topics define research development. A separate issue is the generalization of the R 3 case.