Recovering the polyhedral geometry of fragments

TL;DR

This paper introduces an algorithm to determine the polyhedral geometry of rock fragments from 3D scans, enabling reliable analysis of their shapes and origins.

Contribution

A novel algorithm is introduced for identifying ideal polyhedra in rock fragments using 3D scans, surpassing manual inspection.

Findings

The algorithm successfully identifies face orientations and reconstructs ideal polyhedra from 3D scans.

Benchmarking against human measurements showed the algorithm's effectiveness on 132 fragments.

The method can trace the shape evolution of rounded pebbles back to their original forms.

Abstract

Not only is the geometry of rock fragments often well approximated by ideal convex polyhedra having few faces and vertices, but these numbers carry vital geophysical information on the fragmentation process.

Despite their significance, the identification of the number of faces and vertices of the ideal polyhedron has so far been carried out only through visual inspection. Here, we present an algorithm capable of performing this task in a reliable manner. The input of our algorithm is a 3D scan of the fragment which is a triangulated polyhedron with very large number of faces. Our algorithm performs a systematic simplification using the following steps:•Gaussian smoothing is performed on the spherical histogram of the 3D scans faces to identify the most important face orientations.•Planes carrying the faces of the ideal polyhedron are identified and the ideal polygon is reconstructed•Small faces are removed in a systematic manner

Gaussian smoothing is performed on the spherical histogram of the 3D scans faces to identify the most important face orientations.

Planes carrying the faces of the ideal polyhedron are identified and the ideal polygon is reconstructed

Small faces are removed in a systematic manner

We present two versions of the algorithm that we benchmarked the algorithm against a dataset of human measurements on 132 fragments. Beyond identifying the ideal polyhedral approximation for fragments, our method is also capable of tracing backward the shape evolution of rounded pebbles to their origins.

Image, graphical abstract

Figures

Left, the polyhedron from sample 70 with model (a), right with model (b). The difference is an extra side in algorithm (b) which cuts the wedge short to avoid excess volume. This was not identified as a real polygon side by the expert.Fig 1

Positions of the Fibonacci lattice points on the unit sphere for $[eqn]$ Fig 2The Fibonacci lattice points will serve as the bins of the histogram for the surface polygons. The procedure is the following: For the normal of each surface polygon, we find the bin point which is closest to it (the dot product is maximal) and the value in the bin is increased by the surface area of the polygon (weighted histogram). - 2. Smoothing of the histogram The histogram is very noisy if *N* is large, on the other hand, decreasing *N* gives lower resolution which means that the sides of the reconstructed poly

Confusion matrix between the hand measurement and the algorithm. The numbers indicate the number of pebbles in that category. Left: algorithm (a), right: algorithm (b).Fig 5

Result for object 4 for both algorithms with 7 faces. The face in the front was not identified by the expert.Fig 7 - (3)In some cases, algorithm (a) did not return a fit at all; we list these objects in [Fig. 8](#fig0008). Obviously, these objects are very rounded and although algorithm (b) returns a fit for them their usefulness is questionable. Therefore, we conclude that if algorithm (a) fails to return a polyhedral fit, we may consider the polyhedral approximation of the object is not feasible.Fig. 8Front and side view of objects which fail with algorithm (a), for objects 54, 58, 86 from t

Full text

Specifications tableSubject area:Earth and Planetary SciencesMore specific subject area:FragmentationName of your method:polyhedronizerName and reference of original method:*–*Resource availability:Source code: https://github.com/torokjanos-bme/polyhedronizer,Validation data: https://osf.io/h2ezc/

Background

As observed in [1], the fragmentation of homogeneous, isotropic rock under a uniform stress field (to which we jointly refer as ideal conditions) produces convex polyhedra as fragments and the combinatorial properties (e.g. numbers of faces and vertices) of these polyhedra carry vital (and, sometimes, the only) information about the provenance of the fragments: combinatorial numbers can tell the stress field which produced the fragments. For example, hydrostatic stress splits the rock into a Voronoi-type tessellation with polyhedral fragments with respective average numbers for faces and vertices $[eqn]$ . On the other hand, multiple, successive shear events result in $[eqn]$ . This illustrates that the identification of these numbers is of key importance if we want to use the geometry of the fragments to recover the geophysical process which has produced them. Nevertheless, the identification of the ideal convex polyhedron (carrying these key numbers) is far from trivial, as rock material is neither perfectly homogeneous, nor isotropic, nor is the stress field causing the fragmentation perfectly uniform, once fallen out from the rock fragments are also subject to abrasion processes: in brief, conditions are not ideal nor are the resulting fragments. Our current methods article addresses this problem.

Although the physical fragments we observe resemble convex polyhedra, material imperfections and abrasion often result in partially rounded shapes with indentations and irregularities. Fragments are born from fragmentation processes at time t = 0 and in [1] we have shown that shear and tensile processes create fragment shapes which can be well approximated by convex polyhedra. With time, the shape of the fragment evolves [2], and we assume that our input is the 3D scanned image of a fragment at some t_0_ > 0. We define the ideal polyhedron as the convex polyhedron which serves as the best approximation of the shape at t = 0, neglecting the effects of impurities and later abrasion processes. Recovery of the ideal, convex polyhedral shape can be done by manual count using visual inspection. While such estimates produce good approximations, they, however, depend on the person and thus they are hard to verify and to compare. Since 3D imaging technology is readily available [3], here we present an algorithm and its implementation which recovers the ideal, underlying convex polyhedron in a transparent manner. Our method assigns an ideal convex polyhedron to any scanned particle and provides error estimates to gauge the reliability of the output. We present two versions of our algorithm, both of which reduce the number of faces of the scanned fragments drastically, however, with different emphasis. Version (a) is optimized for finding the correct number of faces of the ideal polyhedron while version (b) is optimized to find the correct volume of the ideal polyhedron. Both versions of the algorithm have two essential control parameters which we tuned, using a large dataset from manual/visual measurements as benchmark. Since the configuration space of convex polyhedral is large (for a polyhedron with E edges this space has d = E-1 dimensions), achieving good match with just two control parameters appears reassuring and suggests that the results are robust and reliable.

Often, scanned fragments are represented by a triangulated surface [3]. In the first step, we construct the convex hull [4]. This has double purpose, first, it ensures that the direction of the surface elements of the object actually align with the surface which is not the case if the fragment have surface irregularities due to material grain structure, second, since we are interested in convex representation of the objects the indentations due to impurities are meaningless for our analysis. The main idea of our method is using statistics of the outer normal of this triangulated convex boundary.

We have validated our method on a dataset [5] produced by manual measurements for [1] (see section Method validation). The dataset contains 132 fragments from the Hármashatárhegy mountain in Budapest, Hungary. The composition of the fragments is dolomite. The fragments were collected from weather induced fractured areas. All fragments were scanned in three dimensions, and we compared the output of our algorithm with the manual/visual data obtained by an expert technician.

Method details

The main purpose of this method is to recover the original (or ideal) convex polyhedral form associated with a scanned object. As noted, the ideal convex polyhedron is directly linked to the generating stress field. We assume that the input is a scanned object S, the surface and the point cloud is associated with a triangulated mesh of its convex hull C(S). The main steps of the algorithm are as follows:

• 1.Spherical histogram: We consider the outer normals of the triangles associated with C(S) and construct a spherical histogram H(C(S)) based on these normals. Additionally*,* we identify a reference point $[eqn]$ , the coordinates of which are the averages of the respective coordinates in our point cloud. The point $[eqn]$ is always inside the convex hull of the object.
• 2.Smoothing of the histogram: We perform Gaussian smoothing [6] with standard deviation $[eqn]$ on the spherical histogram H to obtain the smoothed histogram $[eqn]$ . Note that $[eqn]$ is a scalar function on the sphere.
• 3.Finding the local maxima: We identify the local maxima of $[eqn]$ . We denote the number of the local maxima by the integer-valued function $[eqn]$
• 4.Determining the number of faces: Based on $[eqn]$ , we determine the number $[eqn]$ of the faces of the ideal convex polyhedron. This can be done in two alternative ways:
• (a)We identify $[eqn]$ as the length of the longest plateau of $[eqn]$ or
• (b)We give a priory bound $[eqn]$ and we find the closest value $[eqn]$ of $[eqn]$ Note that this number is subject to change later in the algorithm.
• 5.Determining the normal of the faces of the new polygons: By determining the number $[eqn]$ in this manner, we also determined orientation of the $[eqn]$ faces. Next, we determine the outer normal associated with them using the histogram.
• 6.Determining the distance of the faces from the reference point: We determine the distance $[eqn]$ of each face from the reference point $[eqn]$ .
• 7.The full 3d space is then subdivided by the computed planes, and all closed polyhedra are determined: We construct a hyperplane mosaic [1] $[eqn]$ consisting of the $[eqn]$ selected planes. Note that all cells of $[eqn]$ are convex, some are bounded, some are unbounded. We determine all bounded cells of $[eqn]$ . All bounded cells are convex polyhedra.
• 8.Check the origin: We pick the unique bounded cell which contains the reference point $[eqn]$ . If there is no such cell, then the set of $[eqn]$ planes is not a real representation of a polyhedron, larger value of $[eqn]$ must be chosen.
• 9.Unnecessary faces are removed: In this post-processing step we assume that sufficiently small faces have been created by chipping events where the polyhedron was bisected by a plane. We assign unit volume to the polyhedron with $[eqn]$ faces and under this assumption, for each sufficiently small face we reconstruct the missing part of the polyhedron in the absence of this face. We successively remove faces but only if the volume increases less than a given factor dv.

Our algorithm has two major control parameters. For version (b) we have $[eqn]$ and for both versions we have dv.

One important consequence of our algorithm is that the predicted polyhedra are always simple, namely that exactly three faces and three edges meet at every vertex. This is the result of the reconstruction process which identifies the planes of the sides and not the corners of the polyhedron. In [1] it was shown that this process always leads to a simple polyhedron. A simple polyhedron can be characterized by a single parameter e.g., the F number of faces, since the number V of vertices can be calculated as V = 2F-4. from the simplicity property of the polyhedron and Euler formula.

The difference between the versions (a) and (b) is twofold. First, for our data set algorithm (a) runs 44 min for determining the faces and 20 min to remove excess faces, (b) needed 52 min for the faces and 272 min to remove the unnecessary faces, which is about a factor 5 in running time. Second, the results differed in a way which is illustrated in Fig. 1, where the left figure was created with algorithm (a) the right one with (b). The original pebble has a wedge shape, and on the left-hand side it has a rounded edge due to wear. The examiner did not identify this as a face, therefore his judgment on the values of $[eqn]$ matches with algorithm (a). On the other hand, the volume of model (a) is 35 % larger than the original pebble, whereas for algorithm (b) it is only 20 %. This difference is the result of whether the unrealistic wedge is chipped off or not by the algorithm. The largest volume difference was 660 % in (a), vs. 30 % in (b) for sample 84. We note that it is highly likely that at the birth of the fragment the long wedge was indeed part of the fragment, but it broke after the first few impacts the fragment had.Fig. 1. Left, the polyhedron from sample 70 with model (a), right with model (b). The difference is an extra side in algorithm (b) which cuts the wedge short to avoid excess volume. This was not identified as a real polygon side by the expert.Fig 1

In the following we detail the two versions of the algorithm

Preprocessing We assume that the shape of the three-dimensional body is well approximated by its convex hull. Therefore, our algorithm starts by constructing the convex hull of the original shape. This convex hull is a convex polyhedron defined by a collection of planar faces. In general, the number of faces of the convex hull is orders of magnitude larger than that of the ideal polyhedron. Each face is represented by its vertices and its outward normal vector. From the vertex coordinates the surface area of each face can be computed. These will be the input parameters of our algorithms.

• 1. Spherical histogram The next step is to make a weighted histogram from the normals with the surface area as the weight. The most important problem here is to make reliable binning on the surface of the unit sphere. The Euler angles (θ, φ) are unsuitable for this purpose since they show singular behavior at the north and south pole, and it is difficult to create equal sized bins. Therefore, we used a special tessellation of the surface of the unit sphere which was proven to have 40 % more uniformly distributed cells than any conventional polar tessellation, the Fibonacci lattice [7]. This tessellation uses the angle of the golden ration: $[eqn]$ . Only an odd number of points can be created which we denote as $[eqn]$ . The Euler angles at the meshpoints can be obtained by the following formula (assuming the periodicity of φ).

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ This algorithm gives thus P bin points on the sphere and, according to [7], the distribution of these points on the sphere is close to uniform. Fig. 2. illustrates the positions of the grid points of this algorithm for $[eqn]$ Fig. 2. Positions of the Fibonacci lattice points on the unit sphere for $[eqn]$ Fig 2. The Fibonacci lattice points will serve as the bins of the histogram for the surface polygons. The procedure is the following: For the normal of each surface polygon, we find the bin point which is closest to it (the dot product is maximal) and the value in the bin is increased by the surface area of the polygon (weighted histogram).

• 2. Smoothing of the histogram The histogram is very noisy if N is large, on the other hand, decreasing N gives lower resolution which means that the sides of the reconstructed polyhedra will be very inaccurate in their orientation. Our solution to this problem is to use a smoothing function, or kernel [6]. We chose a Gaussian kernel with the distance between two bin points defined as the shortest path on the surface of the unit sphere. The standard deviation $[eqn]$ of the Gaussian kernel was varied in small steps in the range $[eqn]$ , the lower value resulting in surface distance of a few bin points. A larger parameter range would not add reasonable solutions.We have calculated the histogram for $[eqn]$ in method (a) and for $[eqn]$ in method (b). In all cases we used we increased the standard deviation $[eqn]$ for the Gaussian kernel in 50 equidistant steps from 1/7 to 1. The range of Gaussian kernel standard deviation $[eqn]$ is not considered as a model parameter for two reasons: (a) While we consider the emerging polyhedral structure to be a function of the value of the $[eqn]$ , the final, optimal polyhedron will be selected based on extremal properties of this function (i.e. at the longest plateau). We also mention that the range for $[eqn]$ chosen by our algorithm covers the distance from grid level to radius scale, and therefore all possible scenarios are covered.

Finding local maxima In the Fibonacci lattice we determine the distance $[eqn]$ in such a way that all lattice points have on average 6 others within this distance. All points within $[eqn]$ distance are considered to be neighbors. A bin point is considered a local minimum if all its neighbors contain smaller values.

• 4. Determining the number of faces Method (a): The above histograms and local maximum analysis returns a value $[eqn]$ for the number of faces for all values of N and sigma. We observed that in most cases, if the $[eqn]$ was plotted against $[eqn]$ , the curve displayed a long plateau (see the left panel of Fig. 3.) We can observe that the number of bin points does not have a considerable influence on the results, the plateau is robust and obvious. Despite this we performed the computation for all three N values and identified $[eqn]$ as the value which appeared most in all computations. In some cases, this process does not lead to such an obvious plateau (see Fig. 3 right panel), in such cases Step 9 of the algorithm offers the solution. These cases correspond to more rounded objects, where the number of faces is not obvious, even for a human observer.Fig. 3. Top row: The number of local maxima in the histogram as function of the sigma parameter of the Gaussian kernel. The colors and symbols are for systems with different number of bin points are indicated in the legend in the top left corner of the plot.Bottom row: Red curve with left y axes: The volume of the resulting polyhedron as function of the number of faces in arbitrary units. Green curve with the right y axes: The volume times the number of faces of the polyhedron as function of the number of faces.Left: sample no 9, right: no 8.Fig 3. There were few cases where the faces created in this way did not result in a closed object (especially if one face is rounded and the neighboring faces intersect on the other side of the object unlike the ones in Fig. 1), in these cases we considered the second longest plateau, corresponding to a higher value of $[eqn]$ (see Fig. 4).Fig. 4. The pebbles 9 and 8. The sample on the left is a cuboid, the one on the right-hand side has a more complex structure.Fig 4Method (b): we chose the largest number of detected maxima under 16 for which we found a closed polyhedron.

Determining the normals of the faces of the new polyhedron Once the number F of faces has been determined, we have identified the corresponding bin points on the $[eqn]$ histogram. We identified the normals of the faces of the polyhedron by the vectors corresponding to the bin points.

• 6. Determining the distance of the faces from the reference point To determine the distance of a face from the center of mass of the original object, we must measure the distance of the faces of the convex hull having normal close to the normal of the investigated face. In the specific examples we have considered two unit normals $[eqn]$ to be close if $[eqn]$ (but $[eqn]$ worked almost identically).This step assigns a subset of faces of the original convex hull to each new polyhedral face. For each subset of complex hull faces we computed a weighted histogram, based on their distance to the origin, using the area of the faces as a weight function. We identified the distance by the position of the maximum in the histogram computed in the described manner, using B = 10 equal bins spanning between the largest and the smallest distance. The result is a set of F distances which we denoted by $[eqn]$ .

The full 3d space is then subdivided by the computed planes, and all closed polyhedra are determined In this step we used the algorithm developed in [1] to determine the closed polyhedra. The input of this algorithm is a set of normals $[eqn]$ and the corresponding distances $[eqn]$ of the planes which contain the faces of the polyhedron. It is important that no two normals are parallel (although they can have opposite directions). The algorithm determines the relevant geometric objects in the following order:

• (a)Vertices

• (b)Edges

• (c)Faces

• (d)Polyhedra

• (a) Vertices Any set of three planes meet in one point. We solved $[eqn]$ sets of linear equations and stored all points as possible vertices, together with the identifying indices of the corresponding three planes.

• (b) Edges The set of vertices which were obtained by the intersection of two given planes are collinear and thus they can be ordered using scalar product and two reference points. Once the points are ordered the successive pairs form possible edges. They are stored along with the indices of the two corresponding planes.

• (c) Faces Faces are a series of connected edges forming a closed convex polygon. We start with two adjacent edges that share a common vertex. These two edges determine the plane of the face. Additional edges are added to one end of this polygonal line while ensuring convexity and maintaining the condition that all edges are coplanar: The latter condition is guaranteed by appropriate selection from the list of available edges. On a planar section the edges of the emerging planar mosaic always form a vertex with shape “X,” with exactly four edges that are pairwise aligned. In theory, selecting the edge that maintains convexity with the previous edges is straightforward (by checking the sign of the cross products). However, in computational applications, parallel lines do not exist in an exact sense, therefore we select the edge for which the cross product of the normals of the edges yields the smallest absolute value.This algorithm has two possible stopping scenarios: (i) the next vertex of the newly added edge is the same as the first one of the first edge, then we made a loop, so we found a face, we store it with the plane index and the list of edges. (ii) we do not find any edge which we can connect to the last one in which case this is an open face, and it is discarded.

• (d) Polyhedra Polyhedra created by intersecting planes are necessarily simple polyhedra, meaning that exactly three faces meet at each vertex. We also observe that at every edge, exactly two faces meet. We begin our construction by selecting a vertex and three non-parallel faces that meet at this vertex. Our goal is to complete the polyhedron by successively adding faces to it. It is evident that if a face shares two edges with the incomplete polyhedron, it must be part of the final structure.To facilitate this process, we create a list $[eqn]$ ​, containing the edges of the incomplete polyhedron which belong to a single face. We then search for unused faces that have at least two edges in $[eqn]$ ​. This process continues until either (i) the list $[eqn]$ ​ becomes empty, indicating that a closed polyhedron has been formed, or (ii) all available faces have been considered, but $[eqn]$ ​ remains non-empty. In the latter case, the structure is open to infinity and is therefore discarded.

Check the origin After finding all polyhedra created by the planes determined by the set of $[eqn]$ and $[eqn]$ , we check which one of them contains the origin and return this one. If none of them contains the origin, then we return that this is an invalid set of faces.

• 9. Unnecessary faces are removed The major drawback of using algorithm (a) is that rounded corners and edges may turn into faces which however do not represent any actual face of the original fragment. If the close enough description of the rounded body is a desired outcome one may stop here, however if the interest is in the original, fractured shape from which the scanned object evolved by weathering then step 9 must be performed. In the case of algorithm (b) this step is mandatory.In step 9, all generated faces $[eqn]$ are evaluated for relevance. We check the relevance of $[eqn]$ by tentatively removing $[eqn]$ and computing the access volume $[eqn]$ generated in this step and we call $[eqn]$ relevant if $[eqn]$ is sufficiently large. We implemented this criterion by introducing the threshold $[eqn]$ and for $[eqn]$ we say that the face $[eqn]$ is deemed irrelevant and is removed. The threshold $[eqn]$ ​, which defines the critical volume change, is set to 20 %. This process is repeated until no further faces can be removed.Fig. 3 shows the volume as function of the number of faces of the polyhedron for samples 008 and 009. If the scanned object is close to a polyhedron with small number of faces (sample 009, left panel in Fig. 3) then in this process the volume remains constant until a critical face number when it starts increasing drastically. The resulting F number of faces is very clearly F = 6, both from this analysis and also as concluded from the plateau length in the number of maxima plotted as a function of $[eqn]$ . In the right panel we see few plateaus for object 009. Visual inspection of the object in Fig. 4 also reveals this feature: the object has some small faces which may or may not be included in the final shape.We have found that one can get rid of the parameter $[eqn]$ . If we consider the function *VF (*shown with the green curve in the bottom row of Fig 3.), then the absolute minimum of this function coincides with the longest plateau on the $[eqn]$ histogram. We performed the analysis with this choice for the minimum and its performance is only slightly worse than method (b). It admits mostly a larger number of faces than the other two methods or the hand count, however, this almost always happens in the case of non-polyhedron like objects. This signals that this could be desirable behavior. The supplied code has a switch to use this option.

Summary

Here we summarize the parameters of the algorithms (In brackets we report the values used in this presentation), recommendations for values afterwards:

Algorithm (a):

• •N: the set of number of grid points of the spherical histogram. (values used in example: N = 500,1000,2000), set according to the desired resolution
• •S: The number of values for the sigma of the Gaussian kernel (S = 50), must be set with respect to N to cover grid level to radius with reasonable resolution.
• • $[eqn]$ : The maximum number of sides considered. A reasonable cutoff speeds up the process considerably, for not nice polyhedral objects. ( $[eqn]$ )
• •dn: When making side distance histograms the limit on the 1-cosine difference of the angle between the determined normal and the face of the object (dn = 0.05). Should be adjusted to the resolution and desired number of faces.
• •B: The number of automatic bins for the side distance histogram. (B = 10). May be increased if object’s surface is rough.
• • $[eqn]$ : Relative volume increase allowed by the side removal part ( $[eqn]$ ). This is the ultimate step. which is fast and can be fitted easily. Generally, there is no significant difference between 0.1–0.5.

Algorithm (b):

• •N: The number of grid points of the spherical histogram. (N = 2000), set according to the desired resolution
• •S: The number of values for the sigma of the Gaussian kernel (S = 50), must be set with respect to N to cover grid level and third circle with reasonable resolution.
• • $[eqn]$ : The maximum number of sides considered. ( $[eqn]$ ). Here it is the most important parameter, as an object with this number of sides will be created and then simplified to fit.
• •dn: When making side distance histograms the limit on the 1-cosine difference of the angle between the determined normal and the face of the object (dn = 0.05). Should be adjusted to the resolution and desired number of faces.
• •B: The number of automatic bins for the side distance histogram. (B = 10). May be increased if object’s surface is rough.
• • $[eqn]$ : Relative volume increase allowed by the side removal part ( $[eqn]$ ). This is the ultimate step. which is fast and can be fitted easily. Generally, there is no significant difference between 0.1–0.5.

Method validation

We have validated on a dataset used for the publication [1] and the data can be accessed from https://osf.io/h2ezc/ [5]. It contains 132 fragments from the Hármashatárhegy mountain in Budapest, Hungary. The composition of the fragments is dolomite. The fragments were collected from weather induced fractured areas. All fragments were scanned in three dimension and an expert determined the polyhedral structure.

Fig. 5 displays the confusion matrix of the two algorithms with respect to the hand measurement. The coloring is used to drive the eye. The value of the confusion matrix at position i,j shows the number of objects for which the hand measurement gave i the algorithm j sides. We have also measured two metrics, the average Hamming distance, and the average Euclidean distance. The Hamming distance is 0 if both the hand and the algorithm measurements report the same number of sides, and 1 otherwise, so the smaller the value the better it is. The Euclidean distance is the square root of the sum of the side difference square values also normalized by the number of samples. It is smaller for better matches. We also report the volume fraction, where the volume of the resulting polyhedron is divided by the original one. The results are summarized in the following tables.Algorithm (a)Algorithm (b)Hamming0.570.69Euclidean0.310.09Volume (average ± standard deviation)1.79 ± 0.851.33 ± 0.13Fig. 5Confusion matrix between the hand measurement and the algorithm. The numbers indicate the number of pebbles in that category. Left: algorithm (a), right: algorithm (b).Fig 5

The table shows that although algorithm (a) is better in the Hamming distance metric, namely there are more objects where it matches exactly the experiments, when it fails it is off by a larger amount therefore the Euclidean distance is larger. The biggest difference between the two methods is in the volume of the recovered polyhedra. The biggest misses of algorithm (a) are for objects with rounded edges. In such cases algorithm (b) has superior characteristics.

Limitations

Our algorithm is designed to recover the polyhedral shape of scanned pebbles. It performs as intended if the object does not have rounded faces. Below, we outline some limitations:

• (1)The reconstruction process is not unique, as a small face may result from the original fragmentation or, alternatively, from subsequent weathering which is illustrated in Fig. 6. The algorithm reflects this natural ambiguity; however, it also provides certain parameters to control this, namely its two versions and the parameters $[eqn]$ and $[eqn]$ .
• (2)We have compared the results of the algorithms to human measurements, while we also acknowledge that the latter can also be problematic in some cases. Fig. 7 displays two examples where the two algorithms agree with each other but not with the human expert. Our perspective is that, in these cases, the algorithm outperformed the human expert.Fig. 7. Result for object 4 for both algorithms with 7 faces. The face in the front was not identified by the expert.Fig 7
• (3)In some cases, algorithm (a) did not return a fit at all; we list these objects in Fig. 8. Obviously, these objects are very rounded and although algorithm (b) returns a fit for them their usefulness is questionable. Therefore, we conclude that if algorithm (a) fails to return a polyhedral fit, we may consider the polyhedral approximation of the object is not feasible.Fig. 8. Front and side view of objects which fail with algorithm (a), for objects 54, 58, 86 from top to bottom, respectively.Fig 8
• (4)In summary, the main difference between the two versions of the algorithm is in the handling of small but important faces. From the visual standpoint, merging vertices close to each other could be an option (this is the approach used in the manual reconstruction presented in [1], however, this strategy may produce non-simple polyhedra where the degree of vertices is higher than 3, this it is not compatible with the theoretical framework of the fragmentation process and was therefore discarded. Fig. 6. Results of the two methods for object 40. The left panel shows the results of algorithm (a) with 4 sides, the right with algorithm (b) with 7 faces.Fig 6

Ethics statements

None.

CRediT author statement

János Török: Conceptualization, methodology, code development, writing and editing, Gábor Domokos: Conceptualization, methodology, writing and editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.