{'fraction_non_alphanumeric': 0.03775655800680496, 'fraction_numerical': 0.014048951816485567, 'mean_word_length': 4.308603612351913, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical reconstruction of convex polytopes from directional moments

We reconstruct an n-dimensional convex polytope from the knowledge of its directional moments. The directional moments are related to the projection of the polytope vertices on a particular direction. To extract the vertex coordinates from the moment information we combine established numerical algorithms such as generalized eigenvalue computation and linear interval interpolation. Numerical illustrations are given for the reconstruction of 2-d and 3-d convex polytopes.


Introduction
The reconstruction of the boundary of a shape from its moments is a problem that has only partially been solved. For instance, when the shape is a polygon [23,12], or when it defines a quadrature domain in the complex plane [16], it has been proved that its boundary can exactly be reconstructed from the knowledge of its moments. Both results admit no obvious extension to higher dimension. The reconstruction algorithm for polygons is based on Davis' exact integration formula [6] of a function in the complex plane. For polynomial functions, Davis' formula can be seen as a low dimensional case of identities attributed to Brion [2]. Based on the latter Gravin, Lasserre, Pasechnik and Robins proposed the reconstruction of an ndimensional convex polytope in exact arithmetic [15]. The approach requires finding the roots of a polynomial, whose coefficients are the solution of a linear system, and testing equalities to zero.
Brion's integration formula over a polytope does not relate moment information directly to the vertices of the convex polytope, but rather to the projections of these vertices onto some 1-dimensional subspace. To recover the projections, we recognize an inverse problem that arises in several areas [11,17,19,20,24] and can be solved numerically as a generalized eigenvalue problem.
After recovering the projections of the vertices on various one-dimensional subspaces, remains the problem of matching different projections (in different directions) of the same vertex, with that vertex. In this paper we describe how to solve this issue, without resorting to exact arithmetic. The problem cannot be solved with ordinary interpolation or least squares approximation. But using an interval interpolation technique [26], we understand why we need n + 1 projections (or more) to solve the matching.
Our method is the result of combining techniques from quite different mathematical disciplines: integer lattices, computer algebra, numerical linear algebra, interval methods, inverse problems. The complete algorithm, which we demonstrate in Section 6 and challenge in Section 7, consists of the following steps: 1. The exact number of vertices r is computed from an upper bound R and moments up to order 2R + 1 − n, in a sample of directions.
2. For n + 1 (or more) directions, the projections of the vertices are obtained as the generalized eigenvalues of a structured pair of matrices whose entries are determined from the directional moments up to order 2r − n − 1.
3. Each of these projections is then matched to the corresponding vertex and its coordinates are computed as the coefficients of an n-dimensional interval interpolant.
The different steps in our algorithm involve Hankel matrices, in the singular value decomposition for the computation of r, as well as in the generalized eigenvalue problem delivering the vertex projections. Structured matrices with real elements have condition numbers that grow exponentially with their size [3], and the size of our matrices is determined by the number of vertices of the polytope. In Section 7 we are required to use high precision floating-point arithmetic for the polyhedron with many vertices that represents a brilliant diamond cut.
The paper is organized as follow. In Section 2 we introduce geometric, complex and directional moments together with Davis' and Brion's formulae. In Section 3 we review Prony's method and the related eigenvalue problem to determine the projections of the vertices from the directional moments. In Section 4 we discuss the determination of the number of vertices. In Section 5 we present an algorithm to solve the matching problem. Numerical illustrations are given in the Sections 6 and 7 where we reconstruct 2-d and 3-d convex polytopes.

Geometric and directional moments
In this section we present identities attributed to Brion. These identities are central in [1] for establishing the complexity of the computation of the moments of a polytope. Brion's identities are also at the core of the solution to the inverse problem proposed in [15]. They can actually be seen as a generalisation of Davis' integration formula that was used to solve the shape-frommoment problem in 2D [23,12].
We consider a convex polytope in R n determined by the set of its r vertices V. Abusing the notation, V also denotes the polytope itself.
The geometric moments are The order of the geometric moment m α is |α| = α 1 + . . . + α n . These moments can be expressed as a multivariate polynomial in the coordinates of the vertices. A straightforward formula is given in [27] and a complexity analysis for the computation of those based on Brion's identities is offered in [1] for exact arithmetic.
The moment in the direction δ ∈ R n of order k is where ·, · denotes the usual scalar product in R n .
One can obtain any directional moment of order k from the geometric moments of order k with the multinomial formula Conversely, geometric moments of order k can be obtained from the directional moments of order k in n+k−1 k distinct directions by solving a linear system of equations.
In the context of polygon retrieval (n = 2) from tomographic data in [23,12], geometric moments and then complex moments are computed from directional moments. Complex moments can be understood as moments in the direction δ = (1, i). At the core of this shape-from-moments problem is Davis' integration formula for an analytic function f on the polygon V in the complex where the v ∈ V are here interpreted as complex numbers. Assuming thatv andv are the vertices adjacent to v, the coefficients in Davis' formula are where V v is the oriented area of the parallelogram defined by the vectors with verticesv − v and v − v.
This formula bears a generalisation to any dimension, known as Brion's identities. The formula relates the directional moments with the projections of the vertices. This allows us to work directly with directional moments, which are data that can be deduced from tomographic measurements.
Theorem 2.1. [15] Provided that the orthogonal projections of the r vertices of the convex polytope V on the direction δ are distinct, we have the following equalities where the a v (δ) depend on δ and the adjacent vertices of v in a triangulation of V. Moreover The formula for the coefficients a v (δ) is given in [2, Section 10.3] when V is a simple convex polytope. That is, each vertex in the polytope has exactly n adjacent vertices. Let V v be the set of n adjacent vertices of v. The volume V v of the parallepiped determined by V v is obtained through the determinant of the edge vectors of V v . Then In particular, for a simplex the vertices of which Notice that this is actually a polynomial in δ though we shall use its more compact rational expression.
For a more general convex polytope, one has to consider a partition of the polytope into simplices that does not introduce any additional vertex [2, Theorem 3.1]. The coefficients a v (δ) for the convex polytope is then a sum of its sibblings in the formulae for the simplices. That they do not vanish is proved in [15].
The directions δ ∈ R n to which the theorem applies are those for which u, δ = v, δ for all distinct u, v ∈ V. Those are the generic directions. We examine what happens when δ fails to be generic in this meaning.
On one hand Brion's identities are actually correct for any δ that do not make the denominators of a v vanish. However, if δ is a direction for which the coefficients a v (δ) are well defined but for which there are two distinct vertices u, v ∈ V such that u, δ = v, δ we can write the formula with less than r terms. The linear recurrence introduced in Section 3 is then of order less than r and the associated Hankel matrix is of rank less than r.
In addition, when v 0 , δ = . . . = v p , δ for distinct vertices v 0 , . . . , v p ∈ V that belong to the same simplex of any triangulation, there is a formula similar to (2.1) where the p + 1 terms v i , δ n+k are replaced by terms v 0 , δ n+k , (n + k) v 0 , δ n+k−1 , . . . , (n + k) . . . (n + k − p + 1) v 0 , δ n+k−p . This can be deduced from the rightmost expression in (2.4). The Hankel matrix constructed in Section 3 is still of rank r and admits v 0 , δ as a generalized eigenvalue of multiplicity p + 1. See for instance [9,22].

Recovering the projections of the vertices
In this section we address the problem of retrieving the projections V(δ) = { v, δ | v ∈ V} of the vertices of the convex polytope V from its directional moments m k (δ). We shall recognize an inverse problem that has appeared in several areas, as for instance exponential and sparse interpolation [11,20,24] or Padé approximants [17], beside shape-from-moments [23,12,15]. While [15] approaches the problem with Prony's method, which does not behave so well in floating point arithmetic, we shall favor a formulation in terms of generalized eigenvalues.
The standing assumption is that the projections of the vertices on the direction δ are pairwise distinct. Thus |V(δ)| = |V| = r. Also, we assume in this section that the number of vertices is known. We discuss in next section how this number can be retrieved from only the knowledge of the moments.
From the directional moments (m k (δ)) k we introduce the sequence (µ k (δ)) k∈N of modified directional moments defined by By Theorem 2.1 there exist r non-zero real numbers a v (δ) such that this sequence satisfies The goal is to retrieve the r elements v, δ of V(δ) from (µ k (δ)) k and hence (m k (δ)) k . This is just an instance of the following problem that appears in diverse areas listed above.
Inverse problem: Consider a sequence (µ k ) k∈N such that for some non zero real (or complex) numbers a 1 , . . . , a r and pairwise distinct real (or complex) numbers w 1 , . . . w r , The problem is to find the w i from the knowledge of r and (µ k ) 0≤k≤2r−1 . It can be tackled by Prony's method or as a generalized eigenvalue problem 1 .
First one observes that the sequence (µ k ) k is a solution of a recurrence equation of order r, namely µ k+r = p r−1 µ k+r−1 + . . .
where (−p 0 , . . . , −p r−1 , 1) are the coefficients of the polynomial Applying (3.2) to (µ k ) k∈N for k = 0, . . . , r − 1 leads to the linear system From (µ k ) 0≤k≤2r−1 we can retrieve the characteristic polynomial p of the underlying recurrence by solving the above linear system. The sought numbers w 1 , . . . , w r are the roots of this polynomial. Theses two steps (solving the linear system and computing the roots of the entailed polynomial) is known as Prony's method. It was introduced in [25] and is used in the context of the shape-from-moments problem in [15,23]. The authors of [12,19] introduce a solution in terms of the generalized eigenvalues of a pencil of matrices. It is based on the following facts.
We can recast (3.3) into the matrix equality: The last column of the matrix on the right hand side comes from Equation Let D and W r be respectively the diagonal and the Vandermonde matrices defined by {w 1 , . . . , w r } and appearing in the above equality. The latter can thus be written W r P = DW r . The w i being pairwise distinct, W r is invertible and P W −1 Thus w 1 , . . . , w r are the generalized eigenvalues of the matrix pencils H Computing generalized eigenvalues is a classical problem in numerical linear algebra [13,7,18]. The structured problem we consider here is unfortunately known to be potentially illconditioned. Following [4] we can give an upper bound for the conditioning of the generalized eigenvalue problem as a constant multiplied by the square of the condition number of the Vandermonde matrix W r .
To come back to our initial problem of retrieving V(δ) from (µ k (δ)) k we shall introduce the pencil of Hankel matrices H (1) r (δ) . Its generalized eigenvalues are the elements of V(δ). From those we can construct a matrix of generalized eigenvectors, given by the inverse of the Vandermonde matrix W r (δ). The condition number of W r (δ) is denoted κ(δ). To reduce the conditionning of the generalized eigenvalue problem, which is of order κ(δ) 2 , we shall consider polytopes lying in the unit ball.

Estimating the number of vertices
So far, the number r of vertices has been assumed to be given. But r can also be an unknown of the problem. In this section, we discuss how to numerically retrieve this number from the Hankel matrices H (0) k (δ) formed from the sequence of modified directional moments (µ k (δ)) k in a generic direction δ.
One first observe that H . (4.1) For a generic direction δ, the r elements w i of V(δ) are paiwise distinct and therefore det H Based on this observation, if a strict upper bound R for the number of vertices r is given, then r can be determined as the rank of H (0) R (δ). A caveat is that this matrix may be ill-conditioned. The conditioning of H (0) r (δ) is determined by the conditioning of W r (δ) and A(δ) in (4.1). For this we also refer to the discussion in [11] and [21] that examine the situation in the context of sparse interpolation 2 . The conditioning of the Vandermonde matrix W r (δ) depends on the distribution of the numbers in V(δ) of the projections of the vertices in the direction δ [10]. As for the matrix A(δ), having one of the a i (δ) too small can also lead to an incorrect (numerical) rank for H (0) R . Since we can (even randomly) select multiple directions for the projections, we can retain only those directions for which not both W r (δ) and A(δ) are too ill-conditioned. Alternatively we could apply the rank estimates for Hankel matrices of [5].
Therefore, if we have an overestimation R of the number of vertices we can recover the exact number from the analysis of the numerical rank of H   [13,7]. This is discussed on specific cases in Section 6 and 7.

Reconstruction of the vertices from their projections
In this section we show how to retrieve the set of vertices V from their projections V(δ) onto several directions δ. A difficulty to overcome is that the set of projections we start from are not ordered; we have to match them. This rules out classical interpolation or least square approximation. Our approach is based on the interval interpolation scheme presented in [26]. To simplify our line of arguments on the relative error we assume that the polytope is already contained in the unit ball, a situation to which we shall reduce in practice.
Let ∆ be a set of s > n vectors δ in the unit sphere S n . The assumption is that for each δ ∈ ∆ we have a set of approximationsṼ(δ) for the values of V(δ) = { v, δ | v ∈ V}, the projections of the vertices on the direction δ. The result should be a setṼ of r vectors in R n that consists of approximations for the vertices of the polytope V.
The approximations inṼ(δ) are obtained from the modified directional moments (µ k (δ)) 0≤k≤2r−1 as discussed in Section 3. They are computed as generalized eigenvalues and the conditioning of this problem is given by the square of the condition number κ(δ) of the Vandermonde matrix W r (δ) made of these generalized eigenvalues. Therefore, an element w ofṼ(δ) is actually understood as the center of an interval [w − , w + ] of size 2 (δ) where (δ) should be taken as is the bound on the inaccuracy of the input modified directional moments (μ k (δ)) 0≤k≤2r−1 .
The computation can then continue only when for at least n + 1 directions δ the intervals of V(δ) are disjoint. For simplicity, we assume in the sequel that in all s directions the intervals are disjoint. Note furthermore that with our assumptions, all the elements of V(δ) are in the interval [−1, 1].
For an element inṼ(δ), we do not know which projection is an approximation of which vertex. So we need to find the correct labelling of the projections concurrently with the computation of the vertex coordinates.
An algorithm for the reconstruction of the vertex coordinates can be based on the computation of r linear interval interpolants in n variables, of the form The setṼ of n-tuples of coefficients u used to define those r interpolants are the approximations to the set V of the polytope vertices. The interpolation condition has to reflect the fact that the r functions q u interpolate exactly one w-value perṼ(δ) but s of those values across the sets V(δ). Formally, this can be written as We remark at this point that any set of n projections can be interpolated by a function of the form q u (δ), even if the width of the intervals is zero. At the same time any combination of s projections can be approximated in the least squares sense by a function of the form q u (δ). So none of these classical approaches is very useful in figuring out which projections belong to the same vertex v ∈ V. But an interval interpolant through at least n + 1 disjoint intervals does the job: the nonzero interval width compensates for overdetermining the linear interpolant by at least one interpolation condition. The interval interpolant is stringing the intervals, containing the projections, like beads on the graph of the interpolating function. Any interval, meaning any projection, through which it passes, is marked as belonging to the same vertex.
An interval interpolation problem does not need a priori to have a unique solution: sufficiently small perturbations of the coefficients in the linear form may not violate the interval interpolation conditions. However, a proof is given in [26] for the existence of a unique most robust interval interpolant of the form (5.2). By this we mean the interpolant that stays away as far as possible from violating the interpolation conditions imposed by the interval bounds. The proof is based on the equivalence between the interval interpolation problem and the existence of a Chebychev direction for the polyhedral cone defined by the interval interpolation conditions. A Chebychev direction is a one-dimensional subspace which contains the centers of the inscribed balls with maximal radius (maximal for a certain distance given explicitly in [26]). The existence of a Chebychev direction in our case is equivalent to a convex optimization problem. It follows that the solution is unique.
Because of the labelling problem of the approximate projections, the interval interpolation algorithm becomes a 2-step procedure. In a first step we take the subset∆ of the n best conditioned directions in ∆. Consider the n-tuples in the Cartesian product δ∈∆Ṽ (δ) and select those for which the unique (non-interval) interpolant intersects one interval ofṼ(δ) for each of the s − n remaining directions δ ∈ ∆ \∆. They reflect a correct labelling of the projections. These better conditioned directions∆ are usually near to one another and that drastically cuts down the combinatorial aspect of the procedure by a continuity argument.
Since in this first step the coordinates of the vertices are computed from only n of the s available directions, the obtained values are not maximally accurate. This leads us to the following second step. After ordering the s setsṼ(δ) of r values according to the vertex they are a projection of, the coordinates of each vertex can be computed to maximal accuracy from the total of its s projections. To this end the interval interpolation method can be continued with the computation of the most robust interval interpolants satisfying (5.3). In practice, we obtain it as the least-squares solution to the overdetermined linear system where q u is the linear interpolant defined in (5.2) and w is the corresponding generalized eigenvalues with the correct labelling.

Simulations
We now illustrate the proposed approach for the reconstruction of polytopes from their directional moments. For our simulations we consider centered and scaled polytopes: The origin is the center of gravity of the polytope and the vertices lie in the unit ball. This geometric normalisation corresponds to a transformation on the moments as described in [12].
The vertices of the polytope are to be reconstructed from directional moments. The proposed reconstruction of a convex polytope in dimension n requires directional moments in at least n + 1 generic directions. The order of the moments required in each direction then depends on the number r of vertices: we need directional moments up to order 2r − n − 1. In dimension 2 this can be compared to [23], where the complex moments up to order 2r − 3 are obtained from directional moments up to the same order in 2r − 2 directions.
For a given direction δ, the directional moments m k (δ) are computed in double precision and used to form the entries of the pair of Hankel matrices H  r (δ) in at least n + 1 nearby directions δ as in Section 3. Determine the condition number κ(δ) of W r (δ), the Vandermonde matrix formed with these generalized eigenvalues.
3. Recover the verticesṼ from their approximate projectionsṼ(δ) using the interval interpolation technique described in Section 5 with the error estimate based on κ(δ).
The first two steps are performed using standard numerical linear algebra routines from the NAG library through the Maple interface. In particular, Step 1 makes use of the implementation of the Singular Value Decomposition and Step 2 the QZ-algorithm [28].
Step 3 is implemented in Matlab. Every computations in this Section are performed in double precision.
In step 2, we sample a number of directions and retain those for which the condition number of the Vandermonde matrix W r (δ) is the smallest. The condition number indeed depends on the direction : κ(δ) depends on the distribution of the projected vertices [10]. In particular it increases when the projections of two vertices get closer to one another.
Furthermore, to cut down on the combinatorial complexity of the interpolation scheme in the last step, it makes sense to select directions reasonably close to one another. Selecting a generic reference directionδ with a reasonable condition number κ(δ) and other directions in the neighborhood. Note that here non disjoint intervals for the approximations of the projection on a direction δ would induce a poor condition number κ(δ) of the Vandermonde matrix W r (δ).

Reconstruction of a regular hexagon
The regular hexagon (r = 6) is presented in Figure 1(a) with its symmetry axes and their bisectors. The symmetry axes correspond to the non generic directions. The condition number κ(δ) increases drastically when δ approaches these directions.
The number r of vertices is retrieved as the rank of H (0) R (δ), for R large enough and arbitrary directions δ. Here R = 7 is sufficient to reliably analyse the rank from the singular values. These are plotted in Figure 1(c) for three directions picked at random.
In Figure 1(b) the condition number κ(δ) is plotted for 300 equidistant generic directions. A minimal value is reached for directions bisecting two consecutive axes of symmetry. This leads us to choose the directionδ with angleθ = − 5π 12 . It corresponds to one of the bisectors. We then take 4 nearby directionsθ ± 0.05,θ ± 0.10. For each of the 5 directions we construct the pair of Hankel matrices H 6 (δ) and compute their generalized eigenvalues. For these directions κ(δ) is around 300. We therefore take interval of size 2 × 10 −10 according to (5.1).
We compute the interval interpolants from the 5 sets of projections V(θ + k 0.05), k = −2, −1, 0, 1, 2. The relative error on the computed projections is bounded overall by 4.9 × 10 −12 for the 5 chosen directions. The coordinates of the 6 different vertices of the regular hexagon are recovered as the coefficients of the interpolants graphed in Figure 1(d). The relative error on each computed coordinate compared to its true value is bounded overall by 3.1 × 10 −12 .

Reconstruction of a polygon with 12 vertices
In this second simulation, we consider a centered and scaled 12-gon. It is drawn in Figure 2 After inspecting some directions, we choose the reference directionδ withθ = 0.379521 (arrow in Figure 2 12 (δ) whose entries are obtained from the respective modified directional moments (µ k (δ)) 0≤k≤21 . The condition number κ(δ) of the matrix of generalized eigenvectors W 12 (δ) is around 7 × 10 5 for all 5 directions. The relative error on the computed projections compared to their true values is bounded by 5.0 × 10 −6 .
Note that even the two very close vertices in the top right corner in Figure 2(a) are recovered with the accuracy mentioned above. The distance between them is only of the order of 10 −2 .

Reconstruction of a polyhedron with well-distributed vertices
We first consider the polyhedron with 10 vertices represented in Figure 3(a).
We retrieve the number of vertices of the polyhedron by computing the numerical rank of the Hankel matrix H After inspecting several directions, we select (θ,φ) = (−1.256637, 0.261799) for the reference directionδ and 4 other nearby directions where the condition number κ(δ) is of order 10 4 .δ is indicated by an arrow and a dotted line in Figure 3(a). We take the nearby directions as (θ + ε,φ), (θ,φ + ε), (θ − ε,φ − ε) and (θ + ε,φ − ε) with ε = 0.01. For each direction δ, the pair of matrices H  We compute the 10 linear interval interpolants using intervals of width 2 × 10 −7 for the projections. In Figure 3(c) we show one of those surfaces. The cross denotes the approximate locus of the 5 projections obtained as generalized eigenvalues. The coordinates of the 10 vertices of the polyhedron are obtained as the coefficients of the interpolants. The relative error on this final result is bounded by 5.1 × 10 −6 .  Here we examine a polyhedron with 10 vertices and a triangular face of relatively small size (see the top of Figure 4(a)).
In Figure 4(b) we plot the singular values of H (0) 12 (δ) for 3 random directions δ. The gap between the tenth and the eleventh singular value appears for at least one of the directions.

Diamond
As a challenge we choose an actual brilliant cut of a diamond. It is given as a convex polyhedron with 57 vertices and represented in Figure 5. The stone girdle consists of pairs of vertices very close to one another. The number of vertices and the small distance between the projections of the vertices severely impact the condition number κ(δ) of the Vandermonde matrix W 57 (δ): double precision is no longer enough to retrieve sufficiently accurate values for the projections. We rely on the software floats of Maple to provide the needed number of digits for the computations.
To reliably retrieve the number of vertices, we use a precision of 70 digits. Figure 6  We compute the 57 interpolants in double precision using intervals of width 2 × 10 −4 . The coordinates of the vertices are the coefficients in those interpolants. The relative error is bounded by 7.8 × 10 −5 . We plot the error for all the vertices in Figure 7.