Characterization of structures of particles

A methodology for the characterization of particle structures, especially networks, is developed. This scheme combines 3D image analysis techniques with graph theory tools for the simplification of a structure of thick agglomerates to its skeleton. The connectivity graph of the initial structure is compared with the one of the corresponding skeleton, as a measure of simplification. Examples are used to illustrate the effectiveness of our scheme. Particle structures obtained by Brownian Dynamics simulations are characterized qualitatively and quantitatively. Instead of looking at the characteristics of the structure at the level of the individual particles or neighborhoods of particles, our scheme results in quantitative measures of the network, e.g. the number density of branch-points, the degree of branch-points, and the thickness and the orientation of the branches.


Introduction
Suspensions of particles have been a matter of scientific interest since the beginning of the twentieth century [1,2]. These systems have a lot of applications, e.g. in the food industry [3], cosmetics [4], medicine [5], water treatment [6], paint [7], and ceramics [8]. Their microstructure is interesting both from a scientific and an industrial point of view. The microstructure determines the macroscopic properties of the suspension, especially the mechanical ones [9,10]. In the following, we focus on the characterization of the microstructure of this kind of systems, and especially of network structures.
Many techniques have been used for the characterization of the arrangement of the suspended particles in a suspension, often using only two-particle correlations. The most common is the pair-correlation function, g(r) [11][12][13][14][15][16], and related to that the fractal dimension [17,18]. The radial paircorrelation function is a powerful tool that quantifies the distribution of the inter-particle distances, which, however, give an accurate description only if the structure is homogeneous and spherically symmetric. The cylindrical pair-correlation function in two dimensions is another example [19,20]. To quantify the orientation of the connector vector between any two bonded particles, the second Legendre polynomial, S 2 [20], is used, which quantifies the orientation of the bond direction with respect to an external (imposed) direction.
The information obtained by studying the relative arrangement of only two particles is many times insufficient, so a variety of techniques have been developed concerning the nearest neighbors ( O(10) ) of particles. The Voronoi tessellation is a characterization technique, which is used in various kinds of systems [21][22][23][24][25], particularly, for the characterization of glass forming polymers [26] and in the context of self assembly [27]. For the Voronoi polyhedra, one can quantify their coordination number, the so-called isoperimetric quotient, and the principal moments of inertia [28], which in turn are used for defining measures for prolateness, oblateness, and sphericity/uniformity. The Voronoi tessellation, in combination with its dual Delaunay tessellation [29], is used for determining voids, between atoms or particles [30,31]. Furthermore, adaptive common neighbor analysis in combination with neighbor distance analysis are used mainly for crystalline systems [32], where the Voronoi polyhedra are highly perturbed by slight alternations of the coordinates. Another technique that provides information about the heterogeneity of networks particularly on intermediate length-scales is the Minkowski functionals or quermass integrals [15], based on integral geometry. With this technique, the heterogeneity and the pore-size distribution can be analyzed on different length scales. Another analysis technique is based on the gyration tensor [20,33], calculated for certain clusters of particles, as one can extract information about the size, shape [33], and main orientation [20] of the cluster. Particularly, for analyzing the arrangements of particles in colloidal gels, various other techniques have been employed, namely, the distributions of bond angle [13], angles in bonded triangles [13,32], length of straight paths [34], pore size [35], and the contact distributions [16].
All the above methods require the positions of the particles, although they are not always accessible experimentally. For low-volume fractions, the characterization of colloidalparticle structures can be performed by electron microscopy [36][37][38][39] or by static scattering techniques [40,41], where the latter are directly related to the pair-correlation function, g(r). In denser systems, multiple scattering is observed, which needs to be eliminated or at least limited, e.g., by cross-correlation techniques [42,43], index matching [44], or complementary small-angle neutron scattering [45,46]. The pore-size distribution can be measured in solids, and even under certain conditions in fluids (Cryo-SEM) [14,47]. The most detailed local structural and dynamical information in dense systems has recently been obtained with confocal laser scanning microscopy. This direct visualization technique has been used to study colloidal gels and glasses [48][49][50][51] via time-resolved determination of the three-dimensional positions of the particles. Many characteristics of the particle arrangement can be obtained by the techniques summarized before, however, in networks other features are of interest as well. Networks are observed in colloidal suspensions [13], neural systems [52], polymers [53] and elsewhere [54]. The characterization of colloidal networks could give some insight about the microscopic characteristics that govern the properties of network as a whole, i.e. mechanical properties [9], and transport properties through the network itself [55] and through the suspending medium [56]. For networks, the amount of branch-points (BP), the degree of BP (how many branches emerge form a BP), the thickness of the branches, the existence of percolation, and the number of clusters are key features. These features can be assessed qualitatively by visual inspection ("pattern recognition"), however, it is difficult to perform a quantitative analysis of networks with thick branches and BP in an automated way. To identify the essence of the structure by exploring the characteristics of the neighborhoods of particles (Voronoi) or by cluster analysis turns out to be cumbersome, if not impossible, due to the extensive thickness of the structures (too many particles, with Voronoi cells of small volume and isotropic shape); the essential features of the structure are actually blurred by the multitude of particles. It is proposed in this paper to analyze such thick structures by a combination of an image analysis technique called skeletonization [54,57], and a connectivity graph analysis [58]. Analysis techniques for connectivity graphs have been developed and used extensively in the past [59,60], however, rather than simplifying the connectivity graph of the structure, we apply a simplification on the actual structure. Given the 3D image of a structure, the skeletonization technique makes use of a parallel medial thinning algorithm [61] for the reduction/thinning of the structure to an infinitely thin representation. The characterization of osteocytes has already been performed with this kind of techniques [54,57], although it is also used for other applications, e.g., testing the geometry of porous media obtained by 3D printing [62]. Assessment of the local thickness of particle-structures is possible, defined in terms of the sphere of maximum radius that fits in the structure at that point [63]. The connectivity of the final structure can be simplified using tools of connectivity-graph analysis. The simplification can result to a totally simplified graph, highlighting the essential connectivity between the nodes/ particles. As totally simple graph, we define the graph that has zero global and local clustering coefficients.
In this work, an effort towards the complete quantitative characterization of particle networks is made, especially when thick sub-structures are present. The full detail on the particulate level is going to be translated to the network level, leading to a simplified skeleton where BP and branches can be quantified easily. This paper is organized as follows. First, our scheme for structure simplification is introduced in Sect. 2, discussing in detail also the pre-and post-processing of the structure. After introducing characterization tools for the skeletonized network-structure in Sect. 3, the procedure is validated on some characteristic examples in Sect. 4. The effect of resolution when converting the particle-structure into a 3D binary (voxel) representation is examined in Sect. 5. In Sect. 6 configurations obtained by Brownian Dynamics (BD) simulations [64] with different interaction potentials are analyzed and compared. Conclusions are drawn in Sect. 7.

Method of simplification of thick networks
In this section, we introduce a technique for the reduction of a network of particles to the bare skeleton of the actual structure. For that purpose, we employ an already developed method called skeletonization [54]. In this section, we describe the methodology we use, with all the steps presented in the flowchart in Fig. 1. The procedure includes some pre-processing steps, the creation of a 3D binary image, and the actual skeletonization process, where the actual reduction of the structure happens, see Sect. 2.1. The result of this process is a skeleton (backbone of filled voxels) in the form of a 3D binary image. This skeleton is subjected to refining or simplification, see Sect. 2.2, before the actual quantitative characterization of the network takes place.

Skeletonization
In this section, the 3D image creation and skeletonization of our structures is described. Our purpose is to characterize structures of particles obtained by particle-based simulations. However, the method can be also applied to configurations obtained by experimental techniques [36,54]. In the rest of this section, the steps that are necessary for numerically-obtained structures are discussed, and a distinction is made as far as the experimental configurations is concerned.
Initially, an empty 3D mesh tessellating the simulation box is drawn and a blank (all voxels set equal to zero) binary image is created. For the discretization, the voxel size is chosen such that the radius of the particles R p is commensurate to 10 voxels, i.e., R p = 10 a vox , with a vox being the edge length of the voxel. The volume of a particle is represented by the voxels whose centers are separated from the center of the particle by less than 1.10 R p , see Fig. 2b. Additionally, particles that are bonded (condition for a bond between particles i and j is | ij | ≤ 2.20 R p in our case, where ij is the inter-particle vector) are represented by an additional cylinder of diameter 2.20 R p , between the bonded particles, see Fig. 2b. The reason for choosing the diameter of the particles slightly larger than the actual value is to create sufficiently solid necks between bonded particles. Empty voxels between bonded particles are unwanted, as such empty voxels would result in artifacts, since the thinning process (in the skeletonization) is initiated from the empty voxels next to the surface. All the bond-cylinders are introduced to avoid empty voxels between closely packed particles, e.g., close random packing. The intention is that closely packed particles appear as one body in the skeletonization process, to avoid any empty voxels between the particles. Another reason for choosing this bond-length originates from the interaction potential used to generate the structure: the ratio U exv ∕k B T is lower than 10% at a distance of 2.2R p , where U exv is the excluded volume potential and k B T , the thermal energy. In case of polydisperse systems, the above procedure would have to be amended. For any pair of particles with radii R p,i and R p,j , a bond should be identified if the inter-particle distance satisfies the condition | ij | ≤ 1.10 (R p,i + R p,j ) . And, the cylinder representative of the bond should have a diameter equal to the diameter of the smaller of the two particles.
There is still a chance of finding empty voxels between closely packed particles in the form of cavities, a cavity being a domain of empty voxels that is smaller in extent, in every direction, than the size of the periodic simulation box. Cavities inside the structure create a surface of voxels after the skeletonization. Despite the fact that cavities are of large topological significance [65][66][67], in the systems studied in this paper (i.e., colloidal particles structures) only cavities between particles in close distances are expected. As this is not of interest in this paper, these cavities are filled, see Fig. 2c. To assess the severity of filling the cavities, their number density and volume fraction are quantified, before actually applying the filling. Here, N c is the total number of cavities and N p is the number of primary particles, while V c is the total volume of all the cavities and V is the volume of the box. Normalized in this way, the quantities are defined in a way that they are independent of the size of the box. As already mentioned, the steps necessary for processing experimentally obtained structures are also discussed. The only step that is necessary for experimental results is the last one concerning the filling of the cavities. This step is related to the technique with which the structure has been measured, and an appropriate choice needs to be made with respect to the filling of cavities. The reasons cavities are unwanted during the skeletonization process were mentioned earlier.
Another experimental issue is the noise of the image, and small artefacts created by the method with which the image is obtained [54].
After representing the actual structure as a 3D binary image, we use a 3D medial axis thinning algorithm [61] for reducing it to a one-voxel-thick skeleton. This technique is used to determine network structures in biological systems [54,57]. It was developed for non-periodic systems. However, we intend to use it for configurations that comply with periodic boundary conditions. To avoid issues of periodicity, we use as an input an image that has an edge length of 2L, where L is the box-length or the primary structure. Our primary box is located exactly in the middle of the image, and it has an extra layer of L/2 thickness in every direction which is created by using periodic images of the primary box. The effect of choosing the thickness of the extra layer different from L/2 is not studied further in this paper, because the structures we will be analyzing have their characteristic features on length scales significantly smaller than L/2, and so this choice is not expected to lead to any limitation. The output of the skeletonization is a 3D binary image, where the skeleton is represented by the filled voxels, which are located on a cubic grid. Finally, the actual skeleton is obtained by deleting the extra-layers of thickness L/2 from the primary box.

Simplification of a skeleton
As already mentioned, the skeleton needs post-processing after the skeletonization. Of highest importance in a network is the existence and location of branch-points (BP). To this end, we identify the neighboring voxels by applying a distance criterion for defining bonds, r cut = 1.05 √ 3 a voxel . To locate the neighboring voxels of a skeleton voxel, one should check if some of the 26 neighboring voxels (implied by the grid) also belong to the skeleton. The criterion thus covers the neighbors with a common face, edge, or corner (maximum distance of √ 3 a voxel ). However, there are some issues arising from this, admittedly natural, choice of a Fig. 2 a Schematic representation of the initial structure of particles (blue), in a cubic grid. b Configuration of the particles (blue circles), and binary representation (squares). The blue squares are part of a particle and the red ones part of a bond. c Same configuration as in (b), where the pixels belonging to a cavity have been filled (green squares) distance criterion. The problem is that if two branches emanate from a common point, there is a high probability that the first voxels of the branches heading out from the BP are also neighbors to each other. This means that there will be small triangles of bonded voxels. This is unwanted, as this would eventually result in more than one branch-point on the length scale of one primary particle's radius, i.e. voxel length. Therefore, one must post-process the skeleton to simplifying it in this respect.
Identification of the BPs proceeds along the following steps. First, we determine the voxels that have more than two bonds, and the connectivity between these voxels, see Fig. 3. If there are connected components in this graph, we replace them by one voxel at the position of the center of mass of the voxels that belong to that component. This new voxel inherits the connectivity/bonds from its "parents", where "parents" are the voxels that it is replacing. Then, we reproduce the skeleton with the updated voxels and bonds. The reproduced skeleton may contain triangles of voxels, which relate to a few voxel branches attached to BPs of the skeleton. If there are triangles in the graph, we replace them with one voxel in their center of mass. This new voxel inherits the connectivity/bonds from its "parents". In the case that long bonds are created by the replacement process, we try to shorten these bonds. For this purpose, we add voxels in the middle of the intervoxel separation vector (bond), if r bond > 2 a voxel , so that r bond ≤ √ 3 a voxel . Finally, one does not want several BPs in a volume occupied by a single primary particle, so an additional step is needed for the proper quantification of the BPs, see Sect. 3. In the end, three smoothing steps are applied to the resulting skeleton, where the position of all the voxels that have two bonds is altered according to where i ± 1 denotes the connected voxels, and w is the weight of the smoothing step, in our case w = 1∕3 . If w = 0 , the position of the voxel moves exactly in between the two voxels, and if w = 1 , the voxel stays at its original position. The reason why one should apply some smoothing steps is that the position of the skeleton voxels is dictated by the grid of the binary image used during skeletonization. The smoothing steps help to get rid of unnatural kinks in the skeleton. However, a large number of smoothing steps would result in shrinkage of closed structures, as for example circles, and other unwanted distortions of the structure.

Quantification of simplification
Before we proceed with applying the simplification, we introduce the global and local clustering coefficients for estimating the complexity of the connectivity graph. The global clustering coefficient is defined by the number of triangles divided by the number of triads present in the graph [68,69]. The local clustering coefficient of a node is the number of connections between its neighbors, divided by the amount of potential connections between its neighbors. It becomes unity on a fully connected graph (every vertex is connected to every other vertex) and has typical values ranging from 0.1 to 0.5 in many small-world real-life networks [70].
By following the simplification procedure presented in Sects. 2.1, 2.2, we obtain a skeleton of the simplest possible connectivity, i.e., there are no parallel edges, triangular connections, etc. All the closed triangles, and interconnected 3-bonded voxels have been simplified. The effectiveness of the simplification is investigated with the global clustering coefficient and the bonds per particle/voxel. The global clustering coefficient, and also the local clustering coefficient of the nodes, has a finite value on the initial connectivity graphs of the structures. However, after the application of our simplification scheme these measures (both global and local clustering coefficient) drop to zero, which corresponds to the disappearance of triangles or complex connections, between the nodes. So our final connectivity is far simpler than the initial one. The bonds per voxel is reduced in comparison to the bonds per particle, as a result of the simplification.

Particle-to-voxel transition
The structures targeted for simplification are structures of particles of certain size. However, for turning the structure to a 3D-binary image, a certain resolution, a voxel , has to be chosen. This resolution should be capable of capturing the characteristics of the structure microscopically, and be able to roughly represent the spherical surfaces in terms of general shape (i.e., multiple voxels per surface of one particle).
To this end, the resolution of the binary image should be much finer than the particle-size, which means that multiple voxels in a row are needed to span the diameter of a primary particle. Conclusively, the resulting simplified skeleton has significantly more voxels in absolute numbers than primary particles. That being said, one should note that the skeletonized structure is less complex than the original one, as the reduction of the clustering coefficients shows (see Sect. 2.3).

Characterization of a simplified skeleton
The structure of the totally simplified, and thus easier to characterize, skeleton described in Sect. 2 is quantified in terms of branch-points (BP) (number density, degree), branches (thickness, orientation), existence of percolation, and number of clusters. Details concerning these characteristics are given in the following. The BPs are identified by the number of their bonds to neighboring voxels in the simplified skeleton. Every voxel bearing three or more bonds is identified as a BP. As already mentioned in Sect. 2.2, it is unphysical to have more than one BP in a volume occupied by a single primary particle, however, this artefact can arise due to the skeletonization procedure (during the reduction of the medial surface to the medial axis, the paths created depend on the resolution and the noise of the medial surface) [61]. To avoid such unphysical situations, we proceed as follows. We identify the BPs that have |r IJ | ≤ 2R p , where indices I, J refer to BPs, and we introduce bonds between these BPs. We neglect the rest of the skeleton (voxels that are not BPs), and then the independent clusters of this new connectivity graph are counted as BPs. This step is needed due to the nature of the technique and the structures of particles. Similarly it has been found [61] that when a noisy binary image is skeletonized, first one obtains the medial surface, and then this can be reduced to the medial axis, although the medial axis for different resolution is influenced by the noise on the surface. Therefore, a plane where all the voxels have one of the three coordinates (x, y, z) equal, will be reduced to one branch, however, the same surface with noisy voxel positions is going to create multiple branches. This problem is more prominent when curved planes are present. In this paper, only the medial axis is used, and not the medial surface. After this BP-reduction step, the number density of branch points per volume is given by where N BP is the absolute number of BPs, and N p again is the amount of primary particles. By way of its normalization, n BP is independent of system size. In addition, also the average degree of the BPs can be studied, ⟨d BP ⟩ . The degree of a BP, d BP , is defined as the number of bonded voxels in the simplified skeleton. For the BPs that are replaced during the final step (several BPs inside a primary particle), we use the average degree of the replaced BPs as the degree of the new BP. The difference of the average degree of the replaced BPs with the actual degree of the new BP is insignificant for the structures studied in this paper. These numbers could differ in some specific cases, e.g. a BP with large volume and few branches emerging from it; in this situation, the averaged degree would be larger due to the complicated connections emerging from the skeletonization of the BP, while the degree of the replaced BP would be the number of branches emerging from it.
The characterization of the branches of a network has two components. The first one concerns quantifying the local orientation of the branch, and the second is the estimation of its thickness. The orientation is defined as where i is neither a BP nor directly connected to a BP, and the notation i ± 2 refers to the next-nearest neighbors along the branch. This definition of the orientation is such that is forwardbackward symmetric. Also, it is good to avoid defining orientations based on the directly bonded voxels, since then only 26 discrete directions (13 orientations) would be possible due to the representation on a cubic grid. Choosing voxels that are four bonds apart alleviates this problem substantially. The resulting orientation is an approximation of the tangent of the skeleton at the specific point. For the estimation of the thickness of a branch, we consider the orientation ̂ i of the branch at that specific point on the skeleton i ; in the plane perpendicular to that orientation, we search (starting from the center, i.e., the point on the skeleton) for the edge of the primary image/structure, an idea similar to [71]. Specifically, a grid with voxel size a voxel is created on the plane. The pixels on the plane/grid obtain their values with respect to the nearest voxel belonging to the initial binary image, see Fig. 4. The nearest to i empty voxel is located, and the thickness of the branch d B (i) = 2 d min,i is calculated for every i, where d min,i is the distance between i and the nearest empty pixel. The average of all the values, i.e. for all voxels belonging to a branch, will give the average branch thickness, ⟨d B ⟩ . This quantity will be measured in units of R p .
The orientation vector defined above and used for the determination of the branch thickness can also be used, obviously, to quantify the orientation distribution, namely in terms of the second Legendre polynomial, S 2 , defined as with orientation vectors ̂ of the skeleton and reference orientation ̂ , both of unit length. S 2 is an average measure of the orientation with respect to a reference orientation ̂ . It takes the value −0.5 if the distribution of the orientation vector ̂ contains only vectors perpendicular to ̂ , it is zero if the distribution of the orientation vector is isotropic, and it equals unity for a distribution of orientation which contains only vectors parallel to ̂ . Often, the reference orientation is taken one of the main directions x, y, z.
One could assume that, once the BP-reduction has been applied, also the branches themselves would need to be reduced. However, since the need for BP-reduction originates from noise in the medial surface, the reduction of branches is not necessary, as they already belong to the medial surface. As a consequence, we argue that the results concerning thickness and orientation of the branches will not be altered by reducing the branches. For more details on this point, see Sect. 5.
The existence of percolation in each direction of the box is also investigated. One tries to establish in-box pairing between any two voxels that are close to two opposite faces of the simulation box ("close to" implying that they have a bond crossing that respective face of the box).
The number density of the number of clusters, n cl , is defined by identifying the number of independent components in a connectivity graph of the skeleton. This quantity is reported again with the same normalization as BP number density, Eq. 4.

Validation of the proposed procedure
Some example structures are used to explore the suitability of the technique. The examples are chosen to prove that the procedure outlined above can handle the issues described in Introduction (i.e., the existence of thick branches, thick BPs). Also, the possible existence of cavities is covered in one of the examples. The conclusions of this section serve as validation of the proposed procedure for reducing particlestructures only, and no further measurements or quantitative results about the resulting skeleton are presented in this Section.
For the example structures, an FCC crystal structure is used to create them. Initially, the whole box is filled in with an FCC crystal structure of points. The different structures are then created with a distance criterion from a line or circle. The resulting structures of particles are a thick chain of particles, two hollow structures: a ribbon with a tunnel in the middle and a chain with a cavity in the middle, and finally a multi-structure of a torus and a rounded cylinder.
All the structures presented in this Section are not obtained by simulations of a physical system. They were created as example structures to which our procedure can be Fig. 4 a Example of skeleton structure (red), where specific voxels of skeleton (black) are indicated together with the orientation vector defined by the skeleton at this point (green). b A zoomed-in view of the indicated voxels, with the plane normal to skeleton-orientation vector applied for testing.All our structures have periodic boundary conditions. However, some of the structures are smaller than the box, i.e. they do not span the box in any direction.

Thick chain
A simple example, c.f. Fig. 5, is used to demonstrate the technique. Indeed, the skeletonization technique manages to reduce a thick chain of particles to a thin line of voxels, as expected. One can observe the final skeleton together with the initial image in Fig. 5b.

Hollow structures
A more complicated example of a ribbon with a tunnel, see Fig. 6, is used to further test the technique. One can observe the final skeleton in Fig. 6b. The structure is successfully simplified, and the closed loop is preserved.
An example of a cylinder with a cavity inside, see Fig. 7, is used to illustrate the importance of filling cavities. The resulting skeletonized structure contains a blob of voxels around the (original) cavity in the middle. This structure is unwanted for the characterization of the skeleton, as explained earlier. This problem is solved by filling the cavities that are not connected to the planes of the box, i.e., by using the procedure described in the previous section.

Torus and rounded cylinder
One can see three more complicated examples in Fig. 8. These examples contain a torus structure (thick ring of particles) and a rounded cylinder structure. These structures can be seen as separate structures in Fig. 8a, in contact at four different points in Fig. 8b, and in total contact in Fig. 8c. One can observe the difference in the final skeleton between the different initial images. In Fig. 8a, one can see the two disconnected components of the skeleton. The first is the simplification of the rounded cylinder structure, and the second of the torus structure. The shape of the skeleton thus properly represents the initial image. In Fig. 8b, one can see the same components of the skeleton inter-connected. The connections are located on the points of contact. In Fig. 8c, one can observe only a thin chain of voxels as skeleton. This happens because the two structures (torus and cylinder) are expanded in such a way, that the whole structure is a big agglomerate of voxels with no empty voxels, neither tunnels nor cavities. This example is also a sign that one should be careful with the creation of the initial image. The increase of the radius of the particles should be done with care, since the structure may be altered as in this case. The solution we propose is, as detailed in the previous section, a slight increase of the radius of the particles to avoid small gaps due to the geometry of our primary particles, the introduction of cylinders to the particles we consider as bonded, and the coverage of the cavities of the structure.

Resolution analysis
In this section, we use a structure of an anisotropic network structure, to examine how the results of the structure quantification change upon variation of the number of voxels used per particle radius. In particular, an initial assessment of the region of trust for the resolution is made by computing the accessible volume. Afterwards, in addition to the resolution levels of 2, 5, 10 reported above, also the resolutions of 15, 20 and 25 voxels per particle radius ( R p ) are studied.
The accessible or free volume of the structure is the space that is not covered by a particle and is quantified in an analytic/exact way by using the method developed by Dodd et al. [72], and is then compared with the computed accessible volume (based on empty voxels) from binary images of different resolutions (2-30 voxels/radius) created from the same snapshot (studied in this section), to determine a region of trust for the resolution, see Fig. 9. It is found that beyond a resolution of 10 voxels/radius, the value for the accessible volume of the discretized image converges to a value really close to the one calculated analytically. The variations after a resolution of 10 voxels/radius are small deviations originating from representing spherical particles on a cubic grid. Therefore, a small deviation (of O(0.01%) ) of the numerical from the analytical result is to be expected, and can not be gotten rid of by increasing the resolution further.
The results concerning the filling of the cavities can be found in Fig 10, for the BP (number density, degree) in Fig  11, for the branches (thickness, orientation) in Fig. 12, and for percolation and number of clusters.
One can observe in Fig 10 that higher resolution results in a higher number density of cavities, because smaller gaps between the particles become apparent at higher resolution. However, higher resolution leads to a lower volume fraction of cavities, because the cavities close to the surface of the structure find paths to connect with the outside; the smaller gaps found at higher resolutions are quite small, leading to a decrease in the volume fraction. The result is that the number density is increased, although the volume fraction decreases and convergence is reached after 10 voxels/radius for the volume fraction and after 20 voxels/radius for the number density of cavities. Figure 11 shows that the number density and average degree of the BPs converge as the resolution is increased. The ratio between subsequent resolutions confirms this fact. The increase observed in number density upon increasing the resolution (at low resolutions) occurs because smaller pores (not cavities) in the structure become visible, resulting in more short branches with simpler connectivity (lower degree). The results of Fig. 11 refer to the case when BPreduction has been applied. To illustrate the significance of Fig. 7 Thick chain with a cavity in the interior. a Original arrangement of particles (dark blue). b Corresponding skeleton (dark blue) if cavity-filling is not employed, with superposed original particles (light blue); voxels with more than two bonds in the skeletonised structure are colored green the BP-reduction, we present in Appendix A a regular grid with a known amount of BPs; it is observed that the number density of BPs increases, and does not converge in the absence of BP-reduction, which is unphysical. The results in Fig. 12 show that the average thickness of the branches remains quite constant for resolutions equal to or higher than 5. In this range, the differences are quite small, namely on the order of R p ∕20 , which is comparable to the pixel-edge length. The orientation becomes stronger the higher the resolution, however, the variations being relatively small. This slight variation in orientation can be interpreted by the fact that more branches are present at higher resolution, and the effect of orientation is more pronounced for more detailed skeletons. This issue could be resolved with removing the smaller branches.
We study percolation and number density of clusters for a resolution of 2, 5, 10, 15, 20 and 25 voxels per particle radius. The structure is percolating in x-, y-and z -directions for every resolution, and the number density of clusters is 4.0 × 10 −3 for a resolution of 2, and 3.0 × 10 −3 for the rest of the resolution values. This suggests that the procedure of getting a binary representation of the actual structure with the bonds being represented by cylinders is robust.
In conclusion, one observes that all the presented measures converge when the resolution is increased. Our technique is producing more nodes in the connectivity graphs due to the finer resolution of the voxels compared to the particles, however, this fact is not influencing our measures or the complexity of our graphs in terms of global clustering coefficient. The choice of resolution depends on the structure characteristics one wishes to resolve. If one wishes to fully resolve the cavities in terms of their number density, one should use a resolution of 20 voxels/radius. However, if one wishes to resolve only the network characteristics, i.e. BPs, degree, thickness, orientation, a resolution of 10 voxels/radius is sufficient. The resolution of 10 voxels/radius is used for the rest of this paper.
The image analysis was implemented in MATLAB (MATLAB 2018b, Mathworks), while the network analysis was implemented in python using the Graph tool package [73]. The analysis of a single structure with volume fraction 30% of 1000 particles and resolution of 10 voxels/radius takes less than 12 min on two Intel i5 (i5-8259U) 2.3 GHz CPU and consumes about 1 GB RAM. It is mentioned that not all parts of the code run in parallel.

Application to real structures
The general scheme presented above is now applied to four different structures of different dynamic systems. The structures are thin chains, thick chains, an isotropic network, and an anisotropic network. The thin chains, thick chains, and the anisotropic network are configurations obtained from BD simulations of dielectric particles in a constant and uniform electric field in the z-direction [74], with dipole-dipole interactions [75,76]. The resulting structures are strongly anisotropic in the direction of the external field. In contrast, the isotropic network was obtained again by BD simulations of particles interacting via short-range DLVO-interactions [13]. An important parameter for these structures is the volume fraction of particles, . The different configurations to be analyzed have = 0.05 for thin chains, = 0.10 for thick chains, = 0.20 for the isotropic network, and = 0.30 for the anisotropic network. In general, one tries to establish a relation between the physics included in the simulation and the structures obtained from that simulation. The method presented in this paper bridges this gap, offering quantitative information about the characteristic features of the structure, based on the Cartesian coordinates of the particles. All the structures that are going to be used have periodic-boundary conditions. It should be mentioned that in the connectivity graphs the effect of the periodic-boundary conditions is large. For example, this can be seen for strings that span the box, which will be represented in the connectivity graph as closed loops, due to boundary-crossing bonds between particles on opposite sides of the simulation box (see Fig. 14b).

Cavities
Before presenting the results about the quantification of the actual structures, the presence of (to-be-filled) cavities is assessed in terms of their number density and volume fraction; the results are listed in Table 1.
We expect the thin chains to have the lowest number density and volume fraction of cavities. This is justified by the fact that this structure is the one with the lowest (0.05) and the least local clustering. The thick chains, in contrast, are expected to have more cavities, although still less than the other structures, i.e., the networks. The isotropic network is expected to have more cavities due to the nature of its short-range potential [13,15], and the higher complexity in terms of the form of particle-paths. The volume fraction of cavities is expected to be strongly dependent on the volume fraction of particles, so the isotropic network is expected to have a lower volume-fraction of cavities than the anisotropic network. The anisotropic network is expected to have the largest number density and volume fraction of cavities, the main reason being the volume fraction. In addition, the nature of the interactions (repulsive in direction perpendicular to the external field direction), could promote largevolume cavities. One can see in Table 1, that there is large variation on the measures depending on the structure, the variation being in agreement with the expectation. The thin chains structure has no cavities. The thick-chains structure has the least (non-zero) amount of cavities, and volume fraction, in comparison to the networks. The isotropic network has a one order of magnitude larger number density and volume fraction of cavities. Finally, the anisotropic network shows the highest volume fraction.

Structures
The structures under consideration are presented in Figs. 13, 14, 15, and 16, for the thin chains, thick chains, isotropic network, and anisotropic network, respectively. In these figures, there are four subfigures for all structures, representing the initial structure of the primary particles (a), the final structure of the skeleton (c), and the connectivity of the particles in a 2D projection of the initial structure (b) and of the final skeleton (d). The connectivity graphs are not related to the actual 3D Cartesian positions of the particles/voxels. For all the structure simplification presented in the following, the term "the skeletonization procedure" refers to the entire procedure discussed above, including filling of the cavities and simplification based on the connectivity graph.
One can observe in Fig. 13a, c that the initial structure is preserved, i.e. it is not altered significantly by the skeletonization procedure. This close resemblance of the skeletonized to the initial structure occurs because the structure itself is simple. The only difference is visible in the connectivity graphs, Fig. 13b, d which shows that the structure is indeed simplified, as desired. It shows that even the few complex connections in the initial structure are transformed to simple ones by our procedure. It should be mentioned that the procedure results in more nodes in the connectivity graph of the skeleton than initially present. This happens due to the detailed voxel representation, i.e., since the side length of a voxel is 20 times smaller than the particle diameter, Sect. 2.4. The simplification is verified by the fact that the global clustering coefficient drops from a value of 0.1224 for the initial structure to 0 in the final skeleton. The number of bonds is 1.93 per primary particle; in contrast, the number of bonds is 1.99 on the level of (filled) voxels. This increase is reasonable since the voxel-level gives more emphasis to the branches (many filled voxels are part of a branch) than to the BPs and branch-ends, as compared to the primary-particle level. Therefore, the finer the voxel-representation, the more the number of bonds per (filled) voxel approaches the value 2 (i.e., branch-dominated).
For the structure with thick chains, Fig. 14a, c show that the initial structure is significantly simplified by the skeletonization scheme. The connectivity graphs, Fig. 14b, d, are simplified in terms of bonds per particle/voxel and the local clustering coefficient [68,77]. The connectivity graph is significantly simplified also visually, since the thick rings in Fig. 14b, representative of thick strings, are transformed into single-voxel thick rings in Fig. 14d. The simplification is verified by the fact that the global clustering coefficient drops from a value of 0.3680 for the initial structure to 0 in the final skeleton. The number of bonds is 4.26 per primary particle; in contrast, the number of bonds is 2.03 on the level of (filled) voxels.
Considering the isotropic network, the complexity of the structure in Fig. 15a still allows a reasonable simplification, see Fig. 15c. The connectivity graphs, Fig. 15b, d, are also simplified in terms of bonds per node, and local clustering coefficient [68,77]. However, the simplification cannot be discerned by visual inspection, since the original structure is quite simple already, and only little changes are to be expected from our scheme. The scheme is essentially producing more voxels in the connectivity graph, whereby thin branches are getting longer in the connectivity graph. Also, bonds are crossing and the graph Fig. 13 Structure of thin chains. Initial configuration in real space (a) and its connectivity graph (b); the ratio of bonds per particle is 1.93. Simplified skeleton in real space (red) (c) and its connectivity graph (d); the ratio of bonds per voxel is 1.99. The color in the connectivity graphs indicate different clusters of particles or voxels. In a the size of the particles is reduced by a factor of 2 to increase the visibility of the structure seems rather complicated, although it is fully simplified in terms of global clustering coefficient, which drops from a value of 0.4519 for the initial structure to 0 in the final skeleton. The number of bonds is 1.98 per primary particle; in contrast, the number of bonds is 1.86 on the level of (filled) voxels. The lower ratio of bonds per voxels in the final state is based on the amount of ends present. The skeletonization produces small branches when a structure has small imperfections (see also [61]), resulting to more ends. In the literature, this effect is dealt with a technique called thresholding [54], where the small branches are basically chopped off.
For the anisotropic network, the simplification of the structure is similar to the case of the isotropic network, see Fig. 16a, c. The skeletonization procedure simplifies the connectivity graph, see Fig. 16b, d, in terms of bonds per node, and local clustering coefficient [68,77]. Our procedure is capable of reducing this complex particle structure to its essence, which was the main goal of this paper. The simplification is verified by the fact that the global clustering coefficient drops from a value of 0.3120 for the initial structure to 0 in the final skeleton. The number of bonds is 6.43 per primary particle; in contrast, the number of bonds is 2.10 on the level of (filled) voxels.
The simplified skeletons we obtain from the four structures presented above are further analyzed, to give a more quantitative characterization of these structures.

Branch-points
The results concerning the number density and degree of BP are reported in Table 2. One can observe the differences between the structures. The thin chains have the lowest number density of BP and the lowest degree. This was to be expected as the structure is without significant bifurcations, see Fig. 13. The BPs that do appear are leading to small branches and no connection paths are present between the Fig. 14 Structure of thick chains. Initial configuration in real space (a) and its connectivity graph (b); the ratio of bonds per particle is 4.26. Simplified skeleton in real space (red) (c) and its connectivity graph (d); the ratio of bonds per voxel is 2.03. The color in the connectivity graphs indicate different clusters of particles or voxels. In a the size of the particles is reduced by a factor of 1.5 to increase the visibility of the structure chains. The degree of the BPs supports this statement, since the value of 3.0 is the absolute minimum value of degree a BP can have (by the definition of a BP).
The structure with thick chains has a number density of BP that is almost twice as large as for the thin chains, although it still has a rather low value with respect to the total number of particles. The average degree is slightly higher than for the thin chains. These observations support the existence of more branches and inter-connections between chains.
The isotropic network has the highest values of all four structures for the number density and for the average degree per BP. In the isotropic network, the number density is one order of magnitude higher than for the chain-structures, implying that the structures have a lot of bifurcations, and thus (on average) shorter branches with their thickness being influenced by the volume fraction.
The anisotropic network has a number density and average degree of BP higher than the thin-and thick-chain structures, but lower than the isotropic network. The anisotropic network has double the amount of branch points as compared to the thick-chain structure, but half as compared to the isotropic network.

Branches
The results concerning the characteristics of the branches, namely their thickness and orientation, are shown in Table 3. The thin chains show the lowest thickness and highest orientation along the z-direction. From the measurements of the thickness we understand that there are mostly one-or two-particle-thick chains. The structure with thick chains shows higher thickness, obviously, meaning there are more two-or more particle thick chains. The difference in thickness between the thin and thick chains is rather low, due to the fact that the thick-chain structure looks more like curved sheets than thick cylinders. The orientation with respect to the z-direction is weaker. The isotropic network shows thicker branches, almost two-particle thick on average. The orientation is indeed random, since all three values S 2,i are close to 0. This was expected as the structures have been Fig. 15 Isotropic network. Initial configuration in real space (a) and its connectivity graph (b); the ratio of bonds per particle is 1.98. Simplified skeleton in real space (red) (c) and its connectivity graph (d); the ratio of bonds per voxel is 1.86. The color in the connectivity graphs indicate different clusters of particles or voxels created with a spherical potential, and there is no preferable direction as for the rest of the structures. The anisotropic network structure shows by far the thickest branches, also because this structure has the highest volume fraction. The orientation is weak with a preference for orientation in the z-direction. Fig. 16 Anisotropic network. Initial configuration in real space (a) and its connectivity graph (b); the ratio of bonds per particle is 6.43. Simplified skeleton in real space (red) (c) and its connectivity graph (d); the ratio of bonds per voxel is 2.10. The color in the connectivity graphs indicate different clusters of particles or voxels. In a the size of the particles is reduced by a factor of 1.2 to increase the visibility of the structure

Percolation and Number of Clusters
The results concerning percolation and the number density of clusters are reported in Table 4. The thin chains percolate only in the z-direction, which is the direction of the external electric field, and has the highest number density of clusters, as was to be expected. The thick chains also percolate the box only in the z-direction. This structure has less clusters than the thin-chain structure because more particles make bonds with each other. Both the isotropic and the anisotropic networks percolate in all three directions. The amount of clusters differs significantly. The isotropic network is a result of a diffusion-dominated process. The number of clusters is the lowest of all structures. So all the particles have found their way to a larger structure.

Conclusions
A procedure has been developed for systematic simplification and quantitative characterization of structures formed by particles, with a particular interest in networks with thick branches and branch points. Starting from a structure consisting of thick components, the simplification provides a thin skeleton with minimum connectivity. The procedure can be described in five steps. First, the structure given by the Cartesian coordinates of the particles including bonds is transformed into a discretized binary 3D image in terms of cubic voxels, from which cavities are eliminated. After applying the actual skeletonization process, the resulting skeleton is simplified with respect to the branch points in terms of connectivity reduction and triangle elimination. Finally, the fully simplified skeleton is characterized. In the skeletonization scheme, several numerical parameters can be tuned, which is discussed below. The choices that must be made in the skeletonization procedure concern the pre-and post-processing steps. For the pre-processing, the resolution must be chosen. This choice should be such that it is fine enough to resolve the characteristics of the structure, i.e. accessible volume, but not making the computational cost intractable, because of too many voxels in a the skeleton analysis. Furthermore, one needs to decide whether to eliminate the cavities. One could choose to apply the skeletonization without eliminating the cavities and deal with them in the post-processing steps. The structures investigated in this paper have a low volume fraction of cavities, and we consider them of minor importance for the network structure; in general, one should be critical about the importance of cavities in the structure [66]. In addition, one can decide about including the cylinders representing the bonds. The definition of whether two particles are bonded is based on the fulfillment of a distance criterion that has to be set. Choosing the critical value for the center-to-center distance slightly above the sum of the two particle radii ensures that particles are bonded when represented in the binary 3D image. The critical values can be rationalized on the basis of the interaction potential used to acquire these structures, as explained earlier. In the determination of the thickness of the branches, one should make choices concerning the orientation of the normal vector of the plane, the resolution of the grid on the plane, and the definition of the thickness in this plane. For the orientation, one could choose to use the orientation vector of the bond, instead of the next-nearest neighbors orientation, to determine the plane. The resolution of the grid on the plane can be finer or coarser. And, the thickness can be defined in ways different from our approach, e.g. identification of two empty voxels in opposite directions of the voxel in study, or determining maximum sphere fitting on this point of the image [63].
Standard methods for particle detection require a resolution of approximately 5-20 voxels per diameter of the visualized object, for extracting the particle coordinates [78,79]. If one wishes to apply the thinning algorithm presented in this paper to experimental data, one should use a slightly higher resolution to properly image the characteristics of the structure, since the measures of the network characteristics (e.g. BPs) can be trusted for a resolution that is higher than 10 voxels/radius, or even 20 voxels/radius for the number of cavities. Experimental imaging often results in different resolution in the different spatial directions, the resolution in depth typically being worse than the lateral resolution. In our technique, there is no fundamental reason for the voxels to be cubic, and so parallelepiped voxels of unequal edge lengths can be used in principle. In case the voxels are strongly anisotropic, e.g. much longer in the depth-direction than laterally, one can choose to subdivide the voxel along its depthdirection, thereby generating voxels that are nearly cubic.
Needless to say that the proposed procedure also has some disadvantages. First, choices need to be made in relation to the pre-and post-processing steps, as discussed above. Being able to make choices is in principle not a problem, however, these choices must be conscious and well-thought of. Second, at higher volume fractions ( ≥ 0.40 ), the technique is rather unsuitable, as the unoccupied space is rather low. And third, the application of the technique to a periodic system is cumbersome. One should be mindful about the how thick an extra layer to choose around the primary box, since this is related to the characteristic length scale of the structure.
This skeletonization and analysis procedure was developed for structures obtained via particle-based simulations (e.g., Brownian Dynamics), although it is fully applicable to experimental structures as well. Also in this case, care should be taken with the creation of the image, since additional preprocessing may be needed. Depending on the form of the initial structure, and the way it was obtained, the image may need cavity filling [54]. An advantage or our scheme is its universality, since it can be applied to structures obtained either experimentally or numerically. Being able to use one and the same procedure for analyzing structures of different origin is beneficial for comparing them consistently. ing the code for the accessible volume calculation available to us in Fig. 9. This research forms part of the research program of the Brightlands Materials Center (BMC, www.brightlandsmaterialscenter.com).
Funding KM has received funding by the Brightlands Materials Center (BMC, www.brigh tland smate rials cente r.com).

Compliance with ethical standards
Conflicts of interest The authors declare that they have no conflict of interest.
Source data and code availability The structures are provided as supplementary material, and the code will be uploaded to github, once approval from the sponsor is obtained.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

A Branch-points on a regular grid
A regular grid of 64 intersections is created with 9460 particles. The branches of it are thick FCC structures like the ones presented in Sect. 4, see Fig. 17a. The grid is analysed in terms of BPs for different resolutions, and the results concerning the number density of BP and degree of BP are presented in Table 5. The results about the skeleton before the final reduction (BPs inside the diameter of a particle) are also presented, n BP,before , d BP,before . One can observe that for the BPs before the BP-reduction, at resolution of 20, additional BPs are identified, however, they Table 5 Effect of changing the resolution, in terms of the number of voxels per particle radius, on the number density n BP and average degree ⟨d BP ⟩ of BPs, for the regular-grid structure. The index "before" refers to before the reduction of BPs inside the diameter of a primary particle Resolution n BP ⟨d BP ⟩ n BP,before ⟨d BP,before ⟩ 2 6.7653 × 10 −3 6.0000 6.7653 × 10 −3 6.0000 5 6.7653 × 10 −3 6.0000 6.7653 × 10 −3 6.0000 10 6.7653 × 10 −3 6.0000 6.8710 × 10 −3 6.0313 20 6.7653 × 10 −3 5.7672 7.9281 × 10 −3 5.4533 Fig. 17 a Initial configuration of regular grid with noise (blue). b Final skeleton of regular grid with noise (red). c The BPs are depicted without the rest of the skeleton for resolution of 20. The red color indicates the new BPs obtained from the BP-reduction step, whereas green represents the BPs present before the reduction step are reduced with the elimination of BPs inside a primaryparticle diameter. In Fig. 17, one can observe the initial and final configuration of the regular grid, as well as the BPs before and after the BP-reduction. It is obvious that BPs are present in closer distances in higher resolution.