The secret of gambling with irregular dice: estimating the face statistics of polyhedra

This paper is devoted to the description of a new geometric algorithm to estimate the "face statistics" of an irregularly-shaped polyhedron, i


Introduction
Gambling is as old as human culture.Dice were used for gambling purposes over 5000 years ago.The die is as old an invention as the wheel.When rolling a die, one is interested in the number displayed on the facet facing up after the die has come to a rest.The exact position and orientation of the die do not matter.
Dice are symmetric, which makes the probabilities associated with each individual face trivially equal.However manipulating the probabilities by breaking the symmetry (loaded dice) is also a classical issue in gambling.At the same time there are unusual, inherently asymmetric dice: one of these is the 100-sided zocchihedron [16].The faces of the zocchihedron are of roughly equal size, nevertheless its shape is not a Platonic or a Catalan solid.It has been found that the probabilities associated with its facets are significantly unequal [14].
Beyond gambling, obtaining the face statistics of a polygonal object also has important industrial applications.Devices manufactured in automated assembly lines consist of parts, which are often available in unsorted bulk.To enable further manipulation, the parts must be oriented by feeders.Typically, polygonal parts are dropped to a horizontal surface (e.g. a tray, a vibratory bowl or a conveyor belt) and settle on one of their facets.After coming to rest, they are further processed by various manipulators.Some of these devices change the orientations of parts, whereas others are designed to sort badly oriented ones for recycling.
Traditional methods of designing an efficient feeder for a new part are slow and expensive.The main goals of researchers in manufacturing during the past two decades have been to facilitate feeder design by algorithmic planning and optimization [1,2,3]; and to develop flexible or universal orienting devices [4].
When orienting parts, it is not enough to get all parts to rest on the same facet, since uniform orientation is required by manipulators.Nevertheless many orienting devices can only rotate parts about a vertical axis without changing the face on which they rest.This fact makes the prediction of the landing probabilities associated with individual stable faces of the parts (i.e.face statistics) an important subtask of part feeder design.
The face statistics of an individual shape can be determined by repeated physical experiments or by repeated direct computer simulation of the falling object.Both methods are too slow for certain applications.In such situations, simple estimators of face statistics can be used instead.[6,8,13,5,10,11,7,9] The accuracies of many estimators have been tested by various authors and relatively good performance has been reported in all cases.However, none of them has been tested against more than a few specific objects or a restrictive class of shapes (e.g.square prisms).
Motivated by the lack of reference data, I recently created an extensive dataset of face statistics by the numerical simulation of random polyhedra falling onto a horizontal surface.This was used to evaluate existing estimators and to test and optimize new candidates.This method helped me to a new estimator with average error 6.8% and an improved version of the estimator with average error 3.6%.The second value is roughly half of the average error of the best available estimator [12].In the current paper I focus on the computational algorithm of the first estimator, and briefly sketch the improved version.A more detailed discussion of the physics aspect of the problem, and the simulation methods can be found in [12].

Concepts and notations
This section is devoted to the description of technical concepts used in the paper.
Polyhedra: The surface of any object can be approximated by a polyhedron.A concave object can be replaced by its convex hull, and a non-triangular face can be divided to triangles.Thus, I restrict my attention to convex, triangulated polyhedra.The center of mass of a concave object usually does not coincide with the geometric center of its convex hull.I take this into account by considering specimens with randomly shifted center of mass.
Pose: I want to predict, which face the object is touching the ground in its final configuration.In other words, I am not interested in the exact position and orientation (i.e. the configuration) of the object, the only information I seek is: which point of the object's surface is vertically under its centre of mass O.A set of configurations that share this property are together called a "pose" of the object.This is why face statistics is often termed "pose statistics" in the literature.There are infinitely many poses (one associated with each point of the object's surface), but typically only finitely many of them are suitable for resting on the ground.Objects with rotational symmetry, such as cylinders are exceptional in that they have infinitely many equilibrium poses.Nevertheless these poses can be treated as identical from the point of view of object manipulation.
Instead of associating poses with points of the object surface, one often takes a reference configuration of the object and projects its surface Π centrally to a unit sphere Σ centred at O. The resulting spherical polygon is denoted by Π s .This process yields a one-to-one correspondence between the poses and the points of the sphere.Hence, poses can also be represented by three-dimensional unit vectors (⃗ u).

Centroid solid angles (CSA):
The solid angle of a surface embedded in 3D space from O (centroid solid angle or CSA) is defined as the area of the central projection of the surface to the previously defined sphere Σ.
Representing surface geometry: To represent the shape of an object, I use two functions, which map the 2-sphere to scalars.The "distance function" r and the "support function" R both show how thick the object is in a given direction.For an arbitrary three-dimensional unit vector ⃗ u, r(⃗ u) is the unique positive scalar for which r(⃗ u)⃗ u belongs to the surface of the polyhedron Π; and R(⃗ u) is the only positive scalar for which R(⃗ u)⃗ u belongs to a plane touching Π and perpendicular to ⃗ u (Fig. 1).Both functions are piecewise smooth, with well-defined gradients at smooth points and interval-valued gradients at nonsmooth points.There is an important connection between the support function and the physical behaviour of the object.If the object is in pose ⃗ u and it touches a horizontal support surface form above, its potential energy is R(⃗ u) (times some constant).Hence, an interval-valued gradient of R(⃗ u) containg the 0 vector corresponds to an equilibrium pose and local minima of R to stable resting poses.
Some properties of the r and R functions: It can be shown that the gradients of r(⃗ u) and R(⃗ u) contain 0 at the same points; furthermore r(⃗ u) ≤ R(⃗ u) with equality at equilibrium poses.Both functions are piecewise smooth.The smooth regions are bounded by spherical polyhedra on Σ.For r, it is Π s ; for R it is the central projection of the dual of Π by polar reciprocity to Σ.This is denoted by Π ds .The level curves of R and r are concentric circular arcs within each smooth region.

Motivation for the estimator
My estimator attempts to capture the statistical properties of the complex dynamics of an object dropped onto a horizontal surface.Then, I translate the physical model into the language of geometry.
Let E(t) = U (t) + K(t) denote the total mechanical energy, the potential energy, and the kinetic energy of the object as functions of time t.The reference level of the potential energy is the support surface.The object may not penetrate into the surface, yielding U (t) ≥ R(⃗ u(t)); at the same time, K(t) is non-negative.These inequalities imply R(⃗ u(t)) ≤ E(t) i.e. those poses for which R(⃗ u) exceeds a certain threshold, are not reachable by the object.We assume conventional rolling or dropping of the die, i.e. initial conditions when the initial mechanical energy of the object is enough to make every pose reachable.However its motion is accompanied by energy absorption, and the set of reachable poses shrinks gradually.It is assumed that the energy absorption happens either continuously or through small steps.The function R has s local minima corresponding to the stable poses.Accordingly, there are s − 1 energy levels at which the number of disconnected components of the reachable set increases by one.The splitting events occur at some of the saddle points of the function R(⃗ u).The fragmentation of the reachable set can be described by a "splitting graph" (Fig. 2) in which nodes represent connected components of the reachable set at various energy levels.Each node is labelled by the list of stable poses contained by that component and by the critical value of R at which the component breaks up to smaller pieces.Pairs of directed edges represent the splitting events.
The total energy and the actual pose of an object together determine a node of the graph.Initially, this node is the one labelled by all stable poses.The state of a bouncing object jumps from one node to another due to energy absorption.This process is modelled by a simple, open-loop Markov chain on the splitting graph.
When a component of the reachable set splits to two parts, the actual pose of the object belongs to one of the two parts, which determines the new discrete state of the object.Common sense suggests that the probability of being in a larger component is higher.Motivated by this fact, I assume that the probabilities of transitioning to the two new states are linearly proportional to the centroid solid angles of the two components at the critical energy level.
This rule is used to calculate the likelihood of reaching each discrete state during the motion of the object.Among these, the probabilities associated with the absorbing vertices of the graph (each corresponding to one stable pose) are the predictions of the estimator.
Figure 3: The spherical polyhedron Π ds (thick lines) and the stable manifolds of the saddle points (medium lines) plotted over the contour plot of R associated with the polyhedron of Fig. 1 .

Estimation algorithm
The estimator performs the following steps.

Critical points of R:
The local minima are found by projecting the centre of mass O normally to the plane of each triangular facet of Π.A projection inside the triangle implies a local minimum of R, i.e. a stable pose.The saddle points and the local maxima are found by similar simple calculations.We do not discuss here the possibility of degenerate critical points, since they occur with 0 probability in the case of random polyhedra.One example of the degeneracy is a projection of O to a facet falling exactly to the edge of that facet, corresponding to a marginally stable equilibrium.
Smooth regions of the function R: These are bounded by the edges of Π ds , which can be found in a straightforward way (Fig. 1).
Partitioning Π ds : In order to separate the disconnected components of the reachable set, I construct the stable manifolds of the saddle points of R with respect to the −gradient(R) vector flow (i.e.those two curves per saddle points, which start at the saddle, follow the gradient upstream and arrive to a local maximum of R).Within each face of Π ds , the stable manifolds follow arcs of great circles, which go through the vertex of Π s dual to the examined face of Π ds .Hence they can be followed across one or more facets to their endpoints at one of the local maximum points of R (Fig. 1).
The stable manifolds cut the facets of Π ds to smaller spherical polygons, which generates a "refined" polyhedron Π dsr .The stable manifolds partition the surface of Π dsr .The facets belonging to each partition can be determined easily.Each partition is the basin of attraction of one local minimum of R with respect to the −gradient(R) vector field.Within each partition, the reachable set (for a fixed energy level) is either empty or it is connected and contains exactly one local minimum of R. Parts of the reachable set, which are in different partitions may or may not be connected.Their connectivity is examined in the next step.
Construction of the splitting graph: For low energy level, every stable pose is in a separate component of the reachable set (unless the stable pose is not reachable).As the energy level is increased, individual components grow and fuse at some saddle points of R. To find these points, I consider all saddle points of R in increasing order .Each saddle point lies on an edge of Π dsr .
If the two adjacent facets are in different partitions, then the components of the reachable set in those two partitions become connected at the energy level determined by the value of R at the saddle point.
Transition probabilities: In order to determine the transition probabilities of the splitting graph, I calculate the CSA of the reachable set within each partition for all energy levels corresponding to the saddle points of R. Each component of the reachable set is a combination of circular sectors and triangles, because the smooth parts of R are bounded by the edges of a spherical polyhedron (see Sec. 2), and the level curves of R within each smooth region are circles.There are closed formulas for the CSA of these shapes [15].
Estimated probabilities: As the splitting graph is acyclic, the probability of reaching one of the absorbing vertices is the product of the probabilities along the unique directed path from the initial vertex to the final one.
A polyhedron with n vertices has O(n) edges and facets; and R has O(1) saddle points for typical polyhedra.Hence all steps of the algorithm run in O(n) time.

Accuracy, comparison with other estimators
A set of 1000 random polyhedra has been generated [12], and their falling motion has been simulated 100 times with random initial conditions, using brute force simulation of the Newtonian dynamics with impacts and friction.These numerical experiments provided reference data on the face statistics of these polyhedra.
The performance of an estimator can be measured by the root mean square deviation where N is the total number of stable poses for all object; s j and e j are simulated and estimated probabilities.The results are influenced by the finiteness of the number of trials in the simulation.This can be taken into account by a small correction of σ [12].
Various estimators from the literature have been tested on this set of objects.Their σ values were between 7.0% and 11.2%.The error of the new estimator is σ = 6.8%, which is not significantly better than the previous ones.Nevertheless, a modified version of the new estimator has also been tested, in which the reachable set was bounded by the level curves of r instead of R. Given the close relationship of r and R, the modified estimator works with a very similar algorithm to the one described in Sec. 4. This is not presented in detail.The predictions of the modified estimator are illustrated by Fig. 4. I have no physical Figure 4: Each individual dot in the plots shows the landing probability associated with a a specific stable pose of a random polyhedron based on the average of 100 simulations (vertical coordinate) vs. the prediction of this probability given by the estimator (horizontal coordinate).Points close to the main diagonal are indicators of a better estimator.Nevertheless, deviation from the main diagonal is caused partly by the finite number of trials in the simulation.
argument why such a modification should be beneficial, however the error of the modified estimator is only 4%.By combining one of the earlier estimators with the new one, the error is further reduced to 3.6%.This extension is not presented here.

Conclusion
The efficient design of industrial part feeders requires methods to determine the face statistics of the parts.Finding the statistics by repeated dynamic simulation takes roughly a few minutes on an average PC.There are applications for which faster methods are needed.For example, only a quick estimator can provide the user of computer-aided design software with real-time information about the behaviour of a part being designed; it can also facilitate the automated optimization of part shape through a large number of iterative steps in order to improve its face statistics.For this reason, quick estimators with low error rates are of practical interest.A number of estimation methods have been developed earlier, but their performances have not been compared systematically.Recently, I generated an extensive dataset by simulation, which was used to compare existing estimators and to develop new ones.In the current paper, a new estimator with linear computational complexity has been introduced.The new estimator outperforms all earlier ones.In addition to immediate applications in manufacturing, simple models of the behavior of rigid particles may contribute in the future to an improved understanding of complex issues of granular mechanics, such as loess formation, shielding of river beds or the brazil nut effect.

Figure 1 :
Figure 1: Middle: example of a polyhedron (in its reference configuration) with definitions of the functions r(⃗ u) and R(⃗ u).The plane touches the surface of the object and it is perpendicular to ⃗ u.Resting on the faces marked by 1 and 2 as well as on some other faces on the rear side are stable.Sides: filled contour plots of the r(⃗ u)) (right) and the R(⃗ u) (left) functions associated with the same polyhedron.Light colors indicate high values.Notice that the smooth regions of r correspond to faces of the polyhedron projected centrally to the sphere.Smooth local minimum points within some of these regions correspond to stable resting poses.The smooth regions of R correspond to another spherical polyhedron, which is dual to the original in the classic geometric sense.The stable resting poses are indicated by non-smooth local minima of R.

Figure 2 :
Figure 2: Left: splitting graph of the polyhedron of Fig. 1.The positions of the horizontal parts of edges show the critical energy levels.The polyhedron has 5 stable poses marked by numbers from 1 to 5, and thus four critical energy levels.For illustration, circles indicate the potential energy levels at which a component of the reachable set disappears.Right: the sphere represents the space of poses with level curves of the R(⃗ u) function.Dark shading shows the reachable poses at one of the critical energy levels (dashed line).