Contact orientation distributions for visualisation of granular fabric

Contact orientation distributions are widely used to describe the fabric of granular assemblies and its evolution under load. Different types of visualisations, mainly histograms, are used in the literature to plot these distributions. While there are different ways to create such histograms, however, there is generally no discussion of how the chosen way affects the visualisation of the distribution and its interpretation. We develop in this paper a novel, rigorous framework for discussing contact orientations and their visualisation through histograms. This allows presenting, for the first time, in a unified way several existing visualisations and explaining how they are computed. We identify first some issues in how existing visualisations represent the main features of the contact distribution. We then exploit the framework to introduce new histogram types that avoid these issues and provide a better insight into the granular fabric.


Introduction
Fabric is a generic term that refers to characteristics of the grain-scale microstructure of granular materials [29], for example particle directions, inter-particle contact normal orientations and void space geometry. It is a fundamental concept in understanding, among others, stress induced anisotropy [e.g. 9,13,27,29,30,32,37,39] due to the continuous formation and loss of force-transmitting contacts within the soil skeleton.
A granular assembly will generally exhibit a high number of contacts, making it difficult to draw conclusions when visualising all contacts. We are therefore interested in statistical representations of the distribution of contact vector orientations. Several ways have been proposed in the literature to graphically represent such orientations. We provide here a, necessarily limited, overview of some relevant works.
Oda [29] used a vertical bar histogram to plot the frequency of contact normals, obtained experimentally through thin sections. Later, Oda et al. [30] used polar histograms to show contact orientations in a 2D analogue material; while not explicitly stated, the histogram bars seem to be scaled by bin angle and by the total number of contacts, so as to be comparable to the contact probability density function (PDF). Calvetti et al. [5] used vertical bar histograms to plot the number of contacts in different directions in a 2D material. Yimsiri and Soga [48] consider an analytical PDF with axial symmetry, visualised using a polar plot. Madadi et al. [27] and Bosko and Tordesillas [4] use a polar plot for the PDF for 2D particle simulations, as done by Zhou and Ooi [52] for a simulation of a single layer of spheres. Fu and Dafalias [13] and Marteau and Andrade [28] use polar histograms for contacts in 2D particle assemblies, while Jiang et al. [19] use polar histograms for contact forces in 2D assemblies.
Many authors [e.g. 12,15,25,26,33,34,49] use polar histograms of the angle of the contact vector projected on a plane, in most cases using three orthogonal planes. Zhao et al. [51] use the same approach to plot the projected angle PDF. Khalili et al. [23] use a scaled vertical bar histogram to show the distribution of the magnitude of the vertical component of contact vectors in a discrete-element study. Khan et al. [24] show a polar histogram of the elevation angle for contacts in a 3D assembly. Wiącek et al. [47] similarly 44 Page 2 of 17 show a polar histogram for a 3D simulation, though without explicitly explaining if the quantity binned is the elevation or some other angle.
Wei et al. [40] use a spherical histogram for contact directions while Kawamoto et al. [22] and Zhao and Zhou [50] use a similar plot for inter-particle forces; such plots provide more information than polar histograms but are much more difficult to visually evaluate, and necessarily hide some of the data when in print. Jaquet et al. [18] and Wiebicke et al. [44] use point plots using a Lambert azimuthal equal-area projection, which are then binned to give a heat-map plot. Similar point plots are used in Wiebicke et al. [43] and heatmap plots in Wiebicke et al. [45].
While some of the above visualisations use a contact probability density function, the majority use histograms, which are therefore the focus of the current paper. The interpretation of these visualisations is not always straightforward, and often little explanation is provided on what exactly a histogram shows and why, or on how the histogram is actually calculated. Importantly, there is usually little discussion on how (possibly implicit) choices made when choosing the visualisation type can affect the resulting plot and its interpretation. For this reason we propose a new framework to describe contact orientation distributions and their visualisation using histograms. This framework allows us to examine existing visualisations, identify some issues encountered when using them, and propose new visualisations that can better represent the main features of the contact distribution.
It must be pointed out that any histogram representation will necessary discard a large amount of detailed information on contact orientations. In this paper we focus on the need to have a clear, common understanding of the interpretation of the information being retained and presented, as histograms are widely used in the literature. We do not consider, however, the separate issue of whether a histogram representation is appropriate and/or sufficient to present all fabric information of interest in a specific problem. It is clear that often histograms will need to be complemented by more detailed quantitative analyses, considering for example the fabric tensor evolution.
Note that the term "histogram" is used in the literature with three different meanings: a) binning, a procedure in which a set of values is divided into bins (i.e. intervals) and the number of elements in each bin is counted; b) a graphical representation of the results of binning; and c) the specific graphical representation of binning using adjacent rectangular bars (or circular sectors, for polar histograms). While it is usual (including in the case of contact orientations) to consider together the binning and its graphical representation, we consider here separately the binning procedure, to better understand existing and newly proposed visualisations.
The structure of the paper is as follows: In sect. 2 we define contact vectors and their orientation. Section 3 presents the three sample data-sets used to demonstrate the binning and visualisations used in the paper, specifically in sect. 4 for one-component visualisations and in sect. 5 for two-component visualisations. Finally, the conclusions are presented in sect. 6.

Contact vectors
Two rigid spheres in contact ( Fig. 1) touch at a single contact point, where they share a common tangent contact plane. These concepts can be generalised for the case of deformable and/or non-spherical particles, despite eventual multiple or conforming contacts. We therefore consider that each contact can be fully described by its contact point and contact plane. This is the usual, tacit assumption in the literature.
As an additional major simplification, also widely implied in the literature, we ignore the contact position (except for the discussion in Appendix A). Each contact is therefore represented only by its contact plane, usually described by a unit contact vector , normal to the plane. While a magnitude may be meaningfully assigned to the contact vector, e.g. to indicate an assumed contact area for deformable particles, this is not considered here.
A unit vector is defined by two independent components: e.g. azimuthal angle and polar angle in spherical coordinates, or azimuthal angle and altitude z in cylindrical coordinates (see Fig. 2). The elevation angle, complementary to the polar angle, is also frequently used in the literature; as it leads to the same diagrams as the polar angle, only with different labelling, it is not used in this paper.
A unit vector defines a point on the unit sphere, described using the same two components. We choose here to represent unit vectors as points on the unit sphere, as this greatly simplifies the discussion of concepts like binning. Each contact plane has two contact vectors with the same orientation (parallel) but with opposite sense (pointing in opposite directions). Therefore, we only need to consider the vector orientation, ignoring the vector sense. To do so, we use only one hemisphere of the unit sphere, mapping points from the other hemisphere to their anti-diametrical ones. The choice of hemisphere is important, as it has a major effect on both the binning and its visualisation. To demonstrate this effect, for the coordinate system shown in Fig. 2, we will consider two choices: mapping to the northern hemisphere ( z ≥ 0 ) and mapping to the eastern hemisphere ( y ≥ 0 ). No contact information is lost by mapping to one hemisphere, as anti-diametrical points correspond to the same contact.
Every hemisphere includes a single great circle, e.g. the equator (at z = 0 ) for the northern hemisphere ( z ≥ 0 ). If a pair of contact vectors (i.e. anti-diametrical points) both belong to the great circle, it is necessary to choose which one will be considered. A robust implementation should consider both points, each weighted by 0.5 (i.e. counting each as half a contact); if only few such points exist, randomly choosing one of the two points will however give reasonable results.

Data-sets
Three data-sets are used in this paper to study how to best represent graphically the orientation of the contact normal vectors, and to demonstrate how different visualisations of the same data may lead to different conclusions. The datasets and the software to visualise them as described in this paper are available online [31].

The uniform data-set
The uniform data-set contains 100,000 randomly generated, uniformly distributed unit vectors. This data-set is used as a first benchmark for the different visualisation methods in this paper, as it shows no preferential orientations of the contacts.
Generating uniformly distributed unit vectors is equivalent to picking points uniformly on the unit sphere [42]. Using spherical coordinates with uniform distributions of ∈ [0, 2 ) and ∈ [0, ] does not lead to uniform distribution of points on the unit sphere, resulting instead in a higher density of points near the poles (Fig. 3a). It is instead necessary to use cylindrical coordinates, choosing uniform distributions of ∈ [0, 2 ) and z ∈ [−1, 1] (Fig. 3b).

The HCP data-set
The HCP data-set is a synthetic data-set representing the contacts of 1000 particles taken from a hexagonal close packing. It is used to compare the different visualisations for a corner case where only several discrete orientations exist. The packing is obtained by alternating two layers (A and B) of particles, as shown in Fig. 4. Every particle has 12 contacts, with azimuthal angle and polar angle given by where = arctan √ 2 is the face-vertex-edge angle of the regular tetrahedron.

The DEM data-set
The DEM data-set is taken from a Discrete-Element Method (DEM) simulation of dry zeolite granules subjected to oedometric (1D) compression. Details about the simulation parameters can be found in [20]. A geometry of a rigid cylindrical cell is created (diameter 15mm and approximate specimen height 15mm), particles are randomly generated simultaneously, within a predefined space (same diameter but larger height) and are then left to settle in the cell before being axially compressed. Experimental results were used to validate the material parameters of the DEM model. The narrow grading, number of particles and initial packing density (medium dense with a porosity of 45%) of the specimen are therefore modelled to match the respective values of the experimental set-up [20,21]. Figure 5 shows the DEM simulation timestep from which the data-set used here is taken.
The DEM data-set represents a more realistic contact orientation distribution, which can be encountered in a granular medium. The particles are placed in an axisymmetric cell, 1, 5, 9 on layer A 3, 7, 11 on layer B loaded axially and subjected to gravity in the same direction; they are also spherical and have a very uniform size distribution. As it will be seen in the following, these aspects have a major influence in the resulting contact distribution. Additionally, the DEM data-set does not include the contacts between the particles and the cylindrical cell, as these are only exactly horizontal or vertical and would strongly influence the distribution of the orientations. There are 1944 contacts in the DEM data-set. This will affect visualisations with larger numbers of bins, where small differences in the data will be amplified due to the relatively low number of contacts in each bin. While many DEM simulations of densely packed granular media would involve a much higher number of contacts, it is often useful to extract fabric information for specific regions within the granular material. The choice of this specific DEM data-set allows us to see how a relatively small number of contacts affects the visualisation. An example with a much higher number of contacts is provided in Appendix B.
The DEM and HCP data-sets include both contact vectors for each contact. To simplify the binning code, both contacts are considered and the contact counts are then halved. The three data-sets (uniform, HCP, DEM) are used to demonstrate the different visualisations presented in sects. 4 and 5.  A crucial point made here is that, even for one-component representations, considering the binning process in terms of binning of contact orientations and not of individual components is necessary to fully understand the resulting plot. We use this viewpoint to discuss critically some visualisations proposed in the literature, and propose new visualisations to address the identified issues.

Distribution of azimuthal angles
In a 2D analysis, contact orientation is described by a single component, usually an angle. Polar histograms (also known as rose diagrams) are therefore widely used to observe trends in the contact orientations [e.g. 13,28,52]. It is therefore reasonable to use the same type of diagram to plot the distribution of the azimuthal angle for contact vectors in 3D.
As in the 2D case, we use here bins of equal azimuthal angle (i.e. spanning equal azimuthal angle intervals). This generates bins on the surface of the unit sphere that are spherical lunes (hereinafter simply "lunes") of equal dihedral angle, defined by two half great circles intersecting at the two poles. To consider orientation only, as discussed in sect. 2, it is necessary to map the contact vectors to a single hemisphere. Figure 6 shows the resulting bins: for the northern hemisphere the azimuthal angle is ∈ [0, 2 ) and the bins are half-lunes (right spherical triangles), whereas for the eastern hemisphere the azimuthal angle is ∈ [0, ) and the bins are lunes. Figure 7 shows the results of the data-sets plotted for both northern and eastern hemisphere mappings. Following usual practice, angles in the figures are given in degrees, while angles in the text are given in radians.
As expected, the uniform data-set shows an almost uniform distribution for the azimuthal angle, as no preferred orientation exists. The slight deviations are due to the finite number of contact orientations in the data-set.
The HCP data-set, when mapped to the northern hemisphere, results in 12 equally spaced discrete azimuthal angles with the same number of contacts for each angle. Using 36 bins for the corresponding polar histogram results in all the discrete azimuthal angles falling on the edges of the bins. Due to finite-precision floating-point calculations, this can lead to all contacts with given azimuthal angle being assigned to either the bin preceding the edge or the bin following the edge. In the HCP data-set, minor adjustments (of the order of magnitude of machine precision) were made where necessary to ensure that contacts are always assigned to the bin following the edge. This results in a regular distribution of values in the polar histogram for both the northmapped and the east-mapped orientations. In the latter case, however, all contacts in the x direction are assigned to the first bin (starting at = 0 ) rather than to the last one (ending at = ).
The DEM data-set shows a relatively uniform distribution of azimuthal angle, a result of the macroscopic axial symmetry of the problem being modelled and the relatively narrow grading of the spherical particles. This does not indicate a completely uniform contact fabric in the azimuthal plane, but only one that is uniform with respect to the polar axis (in this example the z axis, which is the axis of the oedometer cell). This is further discussed in Appendix A.
The choice of the hemisphere on which contact orientations are mapped, and specifically its relation to the azimuthal plane, affects the resulting plot. The "north-mapped" plots in the left column of Fig. 7 are not symmetric with respect to the origin. Similar plots presented in the literature are generally symmetric [52], as they consider both contact vectors for every contact point; in such a case, however, the plots become completely equivalent to the "east-mapped" plots like in the right column of Fig. 7. While for the considered data-sets the deviation from point symmetry is limited, in general the north-mapped plot conveys more information on the contact orientations than the east-mapped one. It may however be difficult to describe the physical meaning of this additional information. The north-mapped plot also shows a somewhat higher deviation from the uniform distribution, compared to the east-mapped one, but this is mainly due to it having twice the number of bins and therefore half the number of contacts per bin (54 vs. 108).

Projection on a vertical plane
A common approach in the literature is to project the contact vector on a reference plane and then measure the angle of the projection with respect to a reference axis [12,15,25,26,33,34,49]. This projected angle is then plotted using a polar histogram. In most cases, three orthogonal planes are used, with the resulting three polar histograms giving different views into the contact orientation.  Figure 8 shows the projected-angle plots for the three data-sets, for north-mapped vectors projected on the xz plane. The uniform data-set yields a uniform plot, as expected. The HCP data-set plot correctly shows an equal number of contacts in and out of the xy plane. It is easy however to misread the plot as indicating that 1/6 of the contacts are vertical, i.e. in the z direction, while in reality there are no such contacts. The DEM data-set shows a somewhat uniform distribution of orientations except for an increased number of nearly horizontal contacts. For the coordinate system in Fig. 2, the angle obtained by choosing xy as the reference plane and x as the reference axis is, by definition, the azimuthal angle. The angle obtained can span a range of either or 2 , depending on the relation Fig. 7 Distribution of azimuthal angles for the three data-sets (one per row), mapping orientations on the northern (left column) and eastern (right column) hemispheres. The radial axis shows number of contacts and the circumferential axis shows azimuthal angle between projection plane and hemisphere used for the mapping, as discussed in sect. 4.1. Though both angle ranges appear in the literature, the relevant works do not include the detailed discussion on mapping presented here, which is needed to explain the resulting plot.
Choosing a vertical reference plane (e.g. the yz or xz plane) to project the contact vector does not yield the polar angle, but an azimuthal angle with respect to a different azimuthal plane. This is an important, yet overlooked, aspect: such projected-angle plots must be considered in terms of how meaningful the specific vertical plane is as an azimuthal plane, not in terms of how useful the out-ofplane angle with respect to the horizontal plane is. In the latter case, a plot of the distribution of the polar angle, as described in sect. 4.3, should be used.
The contact orientation distribution for the uniform and DEM data-sets used in this paper is (approximately) axisymmetric with respect to the z axis. This is also true for the common case of triaxial and uniaxial loading conditions, widely used in the literature to study fabric. In all these cases, axisymmetry means that all vertical planes should give (approximately) the same projected-angle distribution.
Yet, projection to a vertical plane is not appropriate in these cases, as shown in Fig. 9. Here, the north-mapped uniform data-set is split into two halves, one with the most horizontal contacts ( ∕3 < ≤ ∕2 ) and one with the most vertical ones ( 0 ≤ ≤ ∕3 ), and the projected angle distribution for each half is plotted separately. In this case, approximately 1/3 of the contacts in the most horizontal half have a projected angle in the more vertical range. This happens because even the most horizontal contacts, when at a small angle from the y axis, will be projected to the zx plane as more vertical (smaller angle from z in the zx plane). While this effect cannot be directly identified in Fig 8b, it is still present, as will be seen by comparing with polar angle plots in sect. 4.3.
Axisymmetric contact orientation distributions should therefore not be visualised by projection on a plane containing the symmetry axis, as it distorts the axisymmetric features of the contact distribution. An appropriate plot of the polar angle distribution should be used instead, as described next.

Equal-polar-angle binning
Given their use to plot the distribution of azimuthal angles, polar histograms have also been used to visualise the distribution of polar (or elevation) angles [1,24]. The simplest approach in this case is to use equal-width bins for the polar angle. The resulting bins on the sphere are spherical zones   Fig. 10 for northmapped contacts.
The first column in Fig. 11 shows the corresponding diagrams for the three data-sets. Like for azimuthal angle plots, the choice of hemisphere is important, with mapping to the eastern hemisphere conveying in this case more information than mapping to the northern one. For the data-sets considered, however, the east-mapped plot is quite symmetric; therefore Fig. 11 uses north-mapped data-sets, with ∈ [0, ∕2] , resulting in more legible diagrams.
The uniform data-set shows how this binning method, whilst mathematically correct, is clearly misleading, as it seems to indicate a preferential orientation on the azimuthal plane ( = ∕2 ). The reason for this is that these polar histograms use bins that subtend equal polar angles, but these do not correspond to equal subtended solid angles (unlike for the azimuthal angle). This is easily seen in Fig. 10, where the northernmost bin covers a substantially smaller area of the hemisphere compared to the southernmost one, and therefore it will include fewer points of a uniform distribution.
By failing to represent a uniform distribution through a uniform plot, this representation makes it difficult to understand how a given distribution deviates from uniformity, therefore hindering the identification of preferential orientations in the data-set. The vertical-plane projection used in the literature does not have this issue, as it uses equal solid-angle bins; it has however other issues for axisymmetric contact distributions, as discussed in sect. 4.2. For this reason, several new representations using the polar angle are presented in the following sections. For ease of comparison, all visualisations of the three data-sets for the polar angle are shown together (as separate columns) in Fig. 11.

Equal-solid-angle binning
To avoid the issues identified in the previous sections, we propose an equal-solid-angle binning of the polar angle.
Here the polar angle is binned so that bins correspond to equal solid angles. This approach adopts the positive aspects of the previous two binning schemes (projection on a vertical plane and equal-polar-angle binning), while avoiding their negative aspects. On the sphere, the bins are equal-area (therefore equal-solid-angle) spherical zones delimited by circles parallel to the azimuthal plane. Figure 12 shows the resulting bins for north-mapped contacts.
In this case the bins span different polar angles, so this binning differs from the one in sect. 4.3.1. For azimuthal angles, however, an equal-solid-angle binning would be identical to the equal-azimuthal-angle binning in sect. 4.1, since equal azimuthal angles subtend equal solid angles.
The solid angle A i subtended between polar angles i−1 and i (with i = 1 … n and i−1 ≤ i ) is the area of the corresponding spherical zone on the unit sphere, given by: Considering n histogram bins, setting 0 = 0 , n = ∕2 , and requiring that all A i are equal, we easily obtain the, unequally-spaced, bin limits for the polar angle Switching from spherical to cylindrical coordinates, a point on the unit sphere has altitude z = cos . Equation (3) therefore becomes It is therefore possible to obtain the equal-solid-angle binning for the polar angles by equal-width binning of the z component of the (unit) orientation vectors.
The resulting polar histograms are shown in the second column of Fig. 11. The uniform data-set clearly corresponds to a (almost) uniform plot. For this reason, this method is preferable to the two methods mentioned previously. The HCP data-set, however, shows a similar trend with the equal-polar-angle binning method, due to the existence of only two discrete values of the polar angle of the contacts. The DEM data-set, finally, shows two preferential orientations (close to 30 • and horizontal contacts), which were not obvious in the equal-polar-angle binning.
One disadvantage of the proposed binning method is that the bin size increases significantly towards the vertical orientation ( = 0 ). This can lead to some information being hidden, for example a concentration of contacts in nearly vertical orientation.
Another potential issue with this visualisation is that it uses bars of different width (i.e. circular sectors of different central angle), which may seem less intuitive than equal-width polar histograms. Equation (4) shows that it is possible to plot the same binning as a horizontal bar histogram for z with equal bar widths. A somewhat similar binning approach is used by Khalili et al. [23, Fig. 7], though (4) z i = 1 − i∕n Fig. 12 Binning of the northern hemisphere using spherical zones of equal solid angle without a detailed explanation of the plot; additionally, a vertical bar histogram is used for visualisation, with the bar values being scaled to be comparable to a probability density function.
Such bar histogram plots may also be difficult to interpret. An alternative approach is therefore presented next, by introducing contact intensity.

Contact intensity
For any binning scheme, it is possible to obtain a uniform polar histogram plot for a uniform contact orientation distribution, by normalising the number of contacts in each bin by the solid angle subtended. We call the resulting quantity contact intensity (number of contacts per solid angle) by analogy with radiant intensity in radiometry (radiant flux per solid angle) [e.g. 36].
The vertical-plane-projection and equal-solid-angle binning schemes described in sects. 4.2 and 4.3.2 use bins that subtend equal solid angles, so plotting the contact intensity just scales the corresponding plots for the number of contacts. For the equal-polar-angle binning described in sect. 4.3.1, however, plotting the contact intensity results in a polar histogram with equal-width bars and with a uniform plot for a uniform distribution of orientations. The resulting plots are shown in the third column of Fig. 11.
As expected, the uniform distribution results in an almost uniform plot. The HCP data-set results in a plot similar to the ones for equal-polar-angle and equal-solid-angle binning, again due to the existence of only two discrete polar angle values. Interestingly, the DEM data-set now shows three preferential orientations instead of the two shown for equal-solid-angle binning, due to the greater resolution of the plot near the pole (i.e. for small polar angles). The DEM data-set appears to have a higher contact intensity for nearvertical contacts, yet it is important to consider that the nearvertical bin of the intensity histogram corresponds to only 28 contacts (over a small solid angle, hence the high value of the contact intensity), while the near-horizontal bin corresponds to 308 contacts. It is expected, therefore, that the intensity of near-vertical contacts will be affected more by small variations in the data-set, making this method more sensitive to such variations.
Note that the contact intensity is equal to the probability density function (PDF) for contact orientations, multiplied by the total number of contacts. Contact intensity histograms can therefore be considered as discrete versions of PDF plots, or PDF plots can be considered as smoothed versions of the contact intensity histograms; while both plots require binning of discrete contact data, the histogram plot has the advantage of showing the binning used. Additionally, the contact intensity distribution contains information on the total number of contacts, which is missing from the PDF; this is important when comparing contact distributions between different particle assemblies.

Contact intensity using quantiles
In general, the number of contacts in each bin varies, possibly leading to issues with resolution like those described in sect. 4.3.3. To avoid this issue of different resolution of the bins depending on the polar angle, we can use contact intensity in conjunction with bins containing equal numbers of contacts (i.e. quantiles), as shown in the fourth column of Fig. 11 which uses 12-quantiles.
To calculate the n b -quantiles (which use n b bins) we first sort in ascending order the values of the polar angles of all contacts (mapped to one hemisphere), obtaining a sorted list of angles j with j = 1 … n c . We then calculate the bin edges as the average between the last polar angle value of one bin and the first value of the next bin, that is where j i = ⌊n c i∕n b ⌋ is the index of the last value of bin i. The first and last edges are not the minimum and maximum actual angle values, but are set to the limits of the possible angle values, in this case e 0 = 0 and e n b = ∕2 . The count of contacts for bin i is easily computed as c i = j i − j i−1 .
As expected, the uniform data-set results into an approximately uniform plot. The resulting plot is essentially the contact intensity plot with equal-solid-angle binning, i.e. a scaled version of the equal-solid-angle plot. Very small differences can be introduced by the specific way in which the bin edges are computed for the quantiles.
There is no HCP data-set plot for this visualisation, as the presence of only two discrete polar angle values leads to a histogram with most bins having zero width. It is therefore not possible to identify those bins on the plot, and of course it is not possible to scale the relevant counts by the (zero) bin solid angle. Visualisation of contact intensity using quantiles does not therefore work for data-sets with few discrete values of contact orientations.
For the DEM data-set, the bin edges obtained are not very different from the ones for the equal-solid-angle binning. The main difference is for the nearly-horizontal bin (the last 12-quantile in terms of polar angle), showing that it has significantly higher intensity than the remaining vectors. Indeed, while the equal-solid-angle, contact intensity, and contact intensity with quantiles plots of the DEM data-set in Fig. 11 are equally valid representations of the underlying data, a comparison between them shows that they highlight different aspects of the data.

Two-component visualization
A two-component visualisation of contact orientations shows the distribution of either number of contacts or contact intensity as a function of two components, e.g. both azimuthal and polar angle. It therefore requires an appropriate binning on a unit hemisphere (the northern hemisphere, in this section). The visualization presented here has already been proposed in the literature [e.g. 22,46], though it has not been widely adopted. It is however presented here in greater detail, and it is explained in terms of the framework already introduced for the one-component visualisations, thus allowing direct comparison of the underlying concepts.

Equal-area subdivision of the hemisphere
There are different reasonable ways to divide the hemisphere into bins. We use equal-area (therefore equal-solid-angle) bins which, as for one-component visualisations, result in uniform binning of a uniform contact orientation distribution. More specifically, we follow the approach by Beckers and Beckers [3], due to its relatively simple implementation. We first divide the hemisphere into n zones using parallels (circles of constant altitude), and then we subdivide each zone i into m i equal bins using equally-spaced meridian arcs. In this way the binning is performed initially by altitude z (or polar angle ) and then by azimuthal angle , but with different edges for in each zone, since each zone has a different number of bins.
We must therefore define the altitude z i of each parallel, with i = 1 … n − 1 and (starting from the north pole) z 0 = 1 and z n = 0 . Each zone i (with i = 1 … n ) from z i to z i−1 has a surface area: and each of the m i equal bins in the zone i has the same area A * , which we require to be independent of i, given by: Equation (7), together with the values z 0 = 1 and z n = 0 , easily yields: which gives the zone boundaries z i for given number of zones n and number of regions for each zone m i .
The trivial case of a single bin per zone ( m i = 1 ) results in the one-component equal-solid-angle binning described in Sect. 4.3.2, with Eq. (8) becoming Eq. (4).
Choosing m i as a linear function of i and substituting into Eq. (8), we get: For a = 2ā and b = −ā we get the simple form: The total number of bins in this case is given by: Figure 13 shows an example of a binned hemisphere. Plotting directly on the hemisphere the binned results means that half the hemisphere results are hidden. Indeed, though it is possible to create a spherical histogram [8,22], the resulting plots can be difficult to interpret. This is especially the case for non-interactive representations of such spherical histograms (e.g. in journal papers, reports, and presentations), in which a large part of the data is obscured and it is also very difficult to extract quantitative data from the plot.

Contact heat maps
To avoid the issues related to visualising results on the hemisphere, we project the hemisphere to the azimuthal plane. We use a Lambert azimuthal equal-area projection [35,41], in (9)  which a point on the unit sphere with cylindrical coordinates √ 1 − z 2 , and z is mapped to a point on the plane with polar coordinates √ 2(1 − z) and . This projects the hemisphere into a circular disc, the individual zones into annuli defined by the radii r i = √ 2(1 − z i ) , and the individual bins on the sphere into planar bins with the same area. For the zones defined by equation (11), the corresponding annuli on the plane have uniformly spaced radii given by: The factor √ 2 can be dropped, as the scale of the resulting plot is arbitrary.
For small data-sets or data-sets with a very strong preferential direction, contact orientations may be directly plotted onto the Lambert projection without having to bin the projection [e.g. 2], but in general binning is required. This can be visualised using a three-dimensional histogram with vertical bars starting from the Lambert projection by extending each bin vertically relative to the number of contacts in each bin, yet a heat-map plot is easier to interpret and present in a publication. The resulting plot is a contact-orientation heat-map using a Lambert projection, or simply contact heat-map. Fig. 14 shows the contact heat-maps for the data-sets considered, using the simple form in Eq. (11) with n = 5 zones and ā = 4 bins in the northernmost zone, for a total of 100 bins. In the contact heat-maps each annulus corresponds to a different polar angle range and each section in an annulus to a different azimuthal angle range. Figure 14a shows that the contact heat-map for the uniform data-set is approximately uniform. This is easier to notice if the colour-map starts from zero (Fig. 14b), in which case the heat-map has almost uniform color. If the colour scale of the heat-map only covers the actual range of number of contacts per bin, then the heat-map shows strong random spatial variations, which however represent relatively small deviations from the average of 1000 contacts per bin.
The heat-map for the HCP data-set clearly identifies the twelve discrete contact orientations present in the data-set. As most bins have a zero count of contacts, starting the colour scale from zero yields the same plot as having the colour scale cover only the range of values present in the plot.
The DEM data-set heat-map shows that the contact distribution is not uniform, however it is difficult to further interpret the plot to extract useful information regarding the main features of this distribution. This is because, as mentioned in sect. 4.1, the DEM data-set corresponds to a problem with macroscopic axial symmetry; the main features of the contact distribution are therefore expected to be independent of the azimuthal angle. The variation of (13) r i = √ 2 i n results across azimuthal angles is therefore expected to be random, with a relatively high magnitude because of the higher number of bins compared to the one-component visualisation, that leads to a low average number of contacts per bin. It is clear that the two-component visualisation, e.g. using a contact heat-map as described in this section, preserves more information on the contact orientation distribution compared to using one or more one-component visualisations. At the same time, the two-component visualisation weakens the visibility of features that depend specifically on one component; in these cases, the onecomponent visualisation can reduce the random variations in the data-set.

Conclusions
Contact orientation histograms are useful to visually describe the fabric of a granular assembly and its evolution under load. For this reason different types of such histograms are widely used in the literature, yet mostly with little or no explanation of what is exactly being plotted and why. The interpretation of these visualisations is therefore not always straightforward.
To provide more clarity, a rigorous framework is developed and presented in this paper for discussing contact orientations and their visualisation through histograms. Original aspects of this framework are a) the introduction of "mapping to a hemisphere" to discuss different ways of calculating orientation and how these affect the resulting plots; b) the clear distinction between histograms as binning and histograms as visualisation; and c) the identification of bins on the hemisphere even for one-component visualisations, to better understand the binning procedure.
This framework allows presenting in a unified way existing one-component visualisations, with a significant result being the identification of issues with equal-polar-angle binning and, in the case of axisymmetric fabric, with the commonly used vertical-plane-projection plots. Additionally, the framework allows the introduction of new histograms that avoid these issues and provide a better insight into the fabric.
More specifically, we have proposed an equal-solid-angle binning for the polar angle, and introduced the concept of contact intensity, resulting in two additional histogram types (equal-polar-angle polar histogram of intensity, and quantile polar histogram of intensity).
We have also presented two-component binning and the resulting contact-heat-map visualisation. Although already used in the literature, we discuss this visualisation type in greater detail, presenting its derivation and how it fits within our general contact-orientation-histogram framework, highlighting also the relation with one-component histograms. especially the case for more complex cases, including for materials with non-uniform particle size distribution.
The framework presented here for contact orientations could be extended to describe histograms for more complex quantities, such as particle directions (for particles for which a direction can be uniquely defined) or inter-particle forces. The latter case is simple when all forces are normal and compressive, so that the force magnitude is used as a weight during binning, but becomes more complicated in the presence of tensile or shear forces.

A An example of position-dependent contact visualisation
In Fig. 15 we propose a new visualisation based on the deviation from the positive radial direction. To calculate this, the contact vectors are projected onto the xy plane, and an angle is calculated from the positive radial direction at the contact point to the projected contact vector. This is the only case in this paper where the contact position is used, as it is needed to determine the positive radial direction. Figure 15 shows that there is a concentration of contacts in the circumferential direction. The nearly uniform distribution of azimuthal angles shown for the DEM data-set in Fig. 7 does not therefore indicate a spatially constant uniform distribution of contact orientations, but only one that is (approximately) axisymmetric. Note that including the contacts with the cylindrical cell would introduce a large number of contacts in the radial direction.

B Evolution of contact orientation
We present here, as an additional example, the use of histograms to track the evolution of the contact orientation distribution for a DEM simulation containing a large number of contacts. The simulation is described in detail by [16] (specifically it is the third of six simulations presented in Table 1 of that paper). The simulation consists in a triaxial compression of a cuboidal sample of 101,623 particles with a maximum-to-minimum particle diameter ratio of 20.6, resulting in more than 180,000 contacts (the exact number varies with strain level).

Fig. 15
Distribution of the angle between the positive radial direction and the projection of the contact vector on the azimuthal plane for north-mapped DEM data-sets. The circumferential axis shows this angle and the radial axis shows the number of contacts Fig. 16 Evolution of polar angle distribution at different strain levels, plotted using 46 equal-solid-angle (left), equal-polar-angle (centre) and quantile (right) bins Figure 16 shows the evolution of the contact orientation polar angle using three different visualisations. All three cases clearly show the change from a uniform distribution to one with more vertical and fewer horizontal contacts. We note that the three visualisations are almost equivalent in this case, as there is a relatively smooth variation of contact intensity with polar angle, without high concentration of contacts in particular narrow regions. The main difference is in the most vertical contacts, where binning by polar angle gives more bins near the vertical direction. This provides more detail, but much of it is random variation as the bins contain fewer points.
The contact heat-maps in Fig. 17 show that the distribution of contacts is mostly, but not completely, axisymmetric.