John Maxwell Campbell
email: jmaxwellcampbell@gmail.com

NEW SERIES INVOLVING HARMONIC NUMBERS AND SQUARED CENTRAL BINOMIAL COEFFICIENTS

Keywords: 1991 AMS Mathematics subject classification. 33E20, 33B15 Harmonic number, Central binomial coefficient, Infinite series, Symbolic computation, Gamma function

which we prove through a variation of our previous integration method for constructing 1 π series.

ABSTRACT. Recently, there have been a variety of intriguing discoveries regarding the symbolic computation of series containing central binomial coefficients and harmonictype numbers. In this article, we present a vast generalization of the recently-discovered harmonic summation formula

∞ n=1 2n n 2 Hn 32 n = Γ 2 1 4 4 √ π 1 - 4 ln(2) π
through creative applications of an integration method that we had previously introduced and applied to prove new Ramanujan-like formulas for 1 π . We provide explicit closedform expressions for natural variants of the above series that cannot be evaluated by state-of-the-art computer algebra systems, such as the elegant symbolic evaluation introduced in our present paper. We also discuss some related problems concerning binomial series containing alternating harmonic numbers. We also introduce a new class of harmonic summations for Catalan's constant G and 1 π such as the series 1. Introduction. In [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], it was noted that the elegant Ramanujanlike formula

(1.1) ∞ n=1 2n n 2 H n 16 n (2n -1) 2 =
12 -16 ln [START_REF] Borwein | Integer powers of arcsin[END_REF] π may be proven through an application of a differential operator with respect to a parameter involved in a known hypergeometric identity; this technique is also used in [START_REF] Choi | Summation formulas involving binomial coefficients, harmonic numbers, and generalized harmonic numbers[END_REF] to prove (1.1). A similar strategy was also applied to prove an equivalent formulation of the known equation in 2017 in [START_REF] Liu | Generalized harmonic number summation formulae via hypergeometric series and digamma functions[END_REF]. Both (1.1) and (1.2) are special cases of a very useful result introduced in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] on the evaluation of series containing factors of the form 2n n 2 H n for indices n ∈ N, which was proven through the use of a beta-like integral transform. While applying parameter derivatives to classical hypergeometric identities only produces specific results on harmonic summations containing squared central binomial coefficients as a factor in the summand, the integration technique presented in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] may be applied much more generally. The power of the integration strategy considered in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] motivates the exploration of further applications of the integral transform given in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], as well as investigations on applications of analogues and variants of this integral transform.

On this latter note, a new integral transform T ln,arcsin given in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] in 2017 was used to evaluate new series for 1 π containing factors of the form 2n n 2 H 2n in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF], letting H 2n = 1 -1 2 + • • • -1 2n denote the alternating harmonic number of order 2n for n ∈ N. The operator T ln,arcsin , which is analogous to the integral transform introduced in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], is used in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] to prove the elegant formula for Catalan's constant G, which is also proven through the use of the generating function for (C n H 2n : n ∈ N 0 ) that was evaluated in [START_REF] Boyadzhiev | Series with Central Binomial Coefficients, Catalan Numbers, and Harmonic Numbers[END_REF], letting C n denote the n th entry in the Catalan sequence. As discussed 3 in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], it is not in general feasible to use known generating functions for sequences such as 2n n H n : n ∈ N 0 to evaluate series with summands containing expressions of the form 2n n 2 H n . In general, the evaluation of infinite series involving products of entries of the harmonic sequence and expressions such as 2n n 2 is difficult, and often leads to surprising and elegant results in the field of classical analysis.

The problem of determining explicit symbolic evaluations for summations containing entries of harmonic-type sequences and central binomial coefficients is a deep and interesting subject that has been explored through the use of many different kinds of classical analysisbased techniques. In [START_REF] Genčev | Binomial sums involving harmonic numbers[END_REF], a variety of infinite summations involving generalized harmonic numbers and central binomial coefficients are evaluated through the use of beta-like integrals. In [START_REF] Chu | Infinite series with harmonic numbers and central binomial coefficients[END_REF], a more abstract way of "depicting" harmonic-like numbers is used, writing

H n = σ 1 1, 1 2 , 1 3 , . . . , 1 n and H 2 n -H (2) n 2 = σ 2 1, 1 2 , 1 3 , . . . , 1 n , letting σ m (x 1 , x 2 , . . . , x n ) = 1≤k1<k2<•••<km≤n x k1 x k2 • • • x km
denote the elementary symmetric function of order m. The authors in [START_REF] Chu | Infinite series with harmonic numbers and central binomial coefficients[END_REF] mainly explore the symbolic evaluation of infinite series involving central binomial coefficients as well as expressions such as σ m 1, 1 3 2 , . . . ,

1 (2n-1) 2
and σ m 1, 1 2 2 , . . . , 1 (n-1) 2 . The results put forth in [START_REF] Chu | Infinite series with harmonic numbers and central binomial coefficients[END_REF] are nicely representative of how mathematical problems concerning the symbolic computation of series involving harmonic-like numbers can be closely connected with seemingly unrelated subjects in the theory of symmetric functions and in number theory.

The series expansions for powers of the inverse sine mapping proven in [START_REF] Borwein | Integer powers of arcsin[END_REF] involve central binomial coefficients and "nested" harmonic-type multisums, including

k-1 n1=1 1 (2n 1 ) 2 n1-1 n2=1 1 (2n 2 ) 2 • • • n N -1 -1 n N =1 1 (2n N ) 2 .
Motivated in part by the main results from [START_REF] Borwein | Integer powers of arcsin[END_REF], and in particular the classical infinite series identity

2 3 arcsin z 2 4 = ∞ n=1 H (2) n-1 n 2 2n n z 2n
known to Ramanujan [START_REF] Berndt | Ramanujan's notebooks. Part I[END_REF], the authors in [START_REF] Mattarei | Congruences for central binomial sums and finite polylogarithms[END_REF] determine congruences for

p p-1 n=1 H n-1 (2) n d 2n n t n (mod p)
for prime numbers p and for special values of d, generalizing congruence results given by Zhi-Wei Sun in [START_REF] Sun | A new series for π 3 and related congruences[END_REF], in which the elegant formula

π 3 48 = ∞ n=1 2 n H (2) n-1 n 2n
n is also introduced. This discussion further illustrates how researching new kinds of subjects concerning summations containing harmoniclike numbers together with central binomial coefficients can lead to unexpected results in both applied analysis and number theory, and surprising connections between these disciplines.

In [START_REF] Guillera | More hypergeometric identities related to Ramanujan-type series[END_REF], new hypergeometric identities related to Ramanujan-like series for 1 π are proven using WZ-pairs, and an equivalent formulation of the beautiful Ramanujan-type [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF] formula

5

was recently proven in [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF] through the use of a parameter derivative applied to a classical 6 F 5 (-1) series identity. These striking results strongly inspire us to explore new techniques for computing series with summands with a factor of the form H n and fixed powers of 2n n in the numerator.

The beautiful formula

(1.4) ∞ n=1 2n n 3 H 2n 2 8n = Γ 6 1 3 (8 √ 3 ln 2 -3π) 24 • 2 2/3 π 4
is proven in [START_REF] Wan | Random walks, elliptic integrals and related constants[END_REF] through the use of special values of the multidimensional integral

(1.5) W n (s) := [0,1] n n k=1 e 2πx k i s dx,
which, as discussed in [START_REF] Wan | Random walks, elliptic integrals and related constants[END_REF], is used in the analysis of uniform planar random walks in the case whereby every step is a unit step, with the direction being random. In particular, the definite integral in (1.5) is equal to the s th moment of the distance in a given random walk, measured from the origin of the plane after a total of n ∈ N steps are taken. The delightful binomial-harmonic series given in (1.4) is proven through an identity for W 3 (s), which shows how series with binomial powers and entries in harmonic-like sequences can have direct applications in the theory of random walks, further motivating the exploration of new applications of the main techniques introduced in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF].

The conjectural equality

(1.6) ∞ n=1 2n n 2 n16 n H 2n = 2 3 ∞ n=0 2n n 2 H 2n (2n + 1)16 n
was introduced in 2014 by Sun, with this conjecture being given alongside a number of conjectural p-adic congruences for truncated sums involving squared central binomial coefficients and harmonic numbers, in [START_REF] Sun | p-adic congruences motivated by series[END_REF]. Considering the striking similarity between the righthand side of (1.6) and the series evaluated in (1.3), it is conceivable that the operator T ln,arcsin and the integral transform from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] could be used to prove the conjecture due to Sun in (1.6) as well as related results in the theory of p-adic congruences. Sun also introduced a plethora of conjectural formulas for infinite sums involving harmonic numbers and centered binomial coefficients in [START_REF] Sun | New series for some special values of L-functions[END_REF]. From the preceding discussion, we are highly motivated to pursue research endeavors on further applications of the integration methods from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] and [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF].

The application of parameter derivatives to hypergeometric identities to prove new results on binomial series containing harmonic numbers was recently discussed in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF]. In [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF], formulas for evaluating harmonic summations of the following forms are proven using this method.

∞ n=1 (2a) n (1 -2a) n (n!) 2 H n 2 n ∞ n=1 (2a) n (2b) n n! a + b + 1 2 n H n n + 1
Letting the parameter a given in the former series be equal to 1 4 , we obtain the series

(1.7) ∞ n=1 2n n 2 H n 32 n
which is evaluated in terms of the gamma function in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF]. Since the integration technique introduced in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] is very useful for evaluating many different kinds of summations involving factors of the form 2n n 2 H n , whereas the summation techniques considered [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF] are given by differentiating specific hypergeometric identities, it is natural to consider the problem of evaluating generalizations of the series in (1.7), through the use of the integration method given in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], which is described in Section 1.1. In particular, for a rational function r(n), it is not obvious as to how to compute generalizations of (1.7) of the form

(1.8) ∞ n=1 2n n 2 H n • r(n) 32 n
following the strategies outlined in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF], or using the generating functions for sequences involving products of harmonic numbers and central binomial coefficients given in [START_REF] Boyadzhiev | Series with Central Binomial Coefficients, Catalan Numbers, and Harmonic Numbers[END_REF].

The infinite series in (1.7) also recently appeared in [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF]. In [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF], it 7 is noted that by Bailey's theorem, we have that

(1.9) 2 F 1 a, 1 -a c 1 2 = Γ c 2 Γ c+1 2 Γ a+c 2 Γ 1-a+c 2 ,
and that (1.7) may be computed in terms of the gamma function by applying the operator ∂ ∂c to both sides of (1.9) in the case whereby a = 1 2 . A formula for a p-adic analogue of (1.7) is also proven in [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF], and many supercongruences for finite sums with central binomial coefficients and harmonic-type numbers are established. The lovely formula

∞ n=1 2n n 3 H n 64 n = 2π(π -3 ln 2) 3Γ 4 3 4
is also proven in [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF] following the "usual" method of applying a parameter derivative to both sides of a known hypergeometric identity; in this case, the classical result known as Dixon's theorem is used. Instead of applying partial derivative operators to classical hypergeometric series identities, in our present article we make use of something of an inverse approach by showing how the integration method explored in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] can be used in a very general way to evaluate series involving harmonic numbers and squared central binomial coefficients.

From the above discussion, we see that the problem of determining closed-form expressions for harmonic summations of the form in (1.8) is an intriguing subject with many potential applications. In our present article, we offer vast generalizations of the formula

(1.10) ∞ n=1 2n n 2 H n 32 n = √ π (π -4 ln 2) 2Γ 2 3 4
that had been noted by Tauraso in 2017 in [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF] and by Nicholson in 2018 in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF] by showing how creative applications of the integration strategy we had previously introduced in [4] can be used to evaluate series of the following forms for z ∈ Z >0 .

n∈N 2n n 2 H n 32 n (n + z) n∈N 2n n 2 H n 32 n (2n -2z + 1) n∈N 2n n 2 H n 32 n (2n -2z + 1) 2
In Section 1.1, we briefly review some preliminary results that we need for the main proofs in our article. In Section 2, we offer a new proof of the evaluation for (1.7) using the integral transform from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], to illustrate the idea of applying this integration technique with respect to summations of the form given in (1.8), and we prove the following new formulas.

∞ n=1 1 32 n 2n n 2 H n n + 1 = 8 - 2Γ 1 4 2 π 3/2 - 4π 3/2 + 16 √ π ln(2) Γ 1 4 2 ∞ n=1 1 32 n 2n n 2 H n 2n -1 = 2Γ 5 4 2 (4 ln(2) -π) π 3/2 + Γ 3 4 (4 ln(2) + π -4) √ 2πΓ 1 4 ∞ n=1 1 32 n 2n n 2 H n (2n -1) 2 = (π -4 ln(2))Γ 1 4 2 8π 3/2 - 2 π (π + 4 ln(2) -6)Γ 3 4 Γ 1 4
.

In Section 3, we show how the proof techniques applied in Section 2 can be generalized, and in Section 4 we introduce a new class of Ramanujan-like series for 1 π using a proof technique from Section 3. Finally, in Section 5, we discuss some future avenues of research related to the main results in our present article.

1.1. Preliminaries. The equality whereby 2 ∞ n=1 (2a) n (2b) n n! a + b + 1 2 n H n 2 n = ∞ n=1 (a) n (b) n n! a + b + 1 2 n H n ,
is proven in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF], for elements a and b in the set of complex numbers such that a + b -1 2 ∈ Z <0 , and this identity is used to prove that the series

∞ n=1 (2a)n(1-2a)n (n!) 2 Hn 2 n is equal to (1.11) √ π 2Γ(1 -a)Γ a + 1 2 ψ (1 -a) + ψ a + 1 2 -ψ (1) -ψ 1 2 , letting ψ(x) = Γ (x) Γ(x)
denote the digamma function, thus leading to the evaluation of (1.7) noted in Theorem 2.1 below. We present a generalization of this theorem in our present article, using the main method applied in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], which may be summarized in Lemma 1.1 below.

As noted above, the elegant formula in (1.10), which has motivated much of the work put forth in our present paper, also appears in [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF] and is proven through a straightforward application of Bailey's theorem. We offer a new proof of (1.10) that is significantly different compared to the proofs of this result from both [START_REF] Tauraso | Supercongruences related to 3 F 2 (1) involving harmonic numbers[END_REF] and [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF]. Variants of our proof of (1.10) may be used to greatly generalize the formula in (1.10) through the use of the following fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], as we later discuss.

Lemma 1.1. Letting f : N 0 → C and g(n) = f (n) 16 n (2n-1) , if f is such that ∞ n=0 (-1) n 1 2 n f (n) x 2n ln 1 -x 2 √ 1 -x 2 is integrable on [0, 1], it follows that ∞ n=0 g(n) 2n n 2 H n is equal to 2 π times the sum of 1 0 ∞ n=0 (-1) n 1 2 n f (n) x 2n ln 1 -x 2 √ 1 -x 2 dx and π ln(2) ∞ n=0 ( 2n n ) 2 f (n) 16 n (1-2n) .
The above lemma often allows us to express an infinite series of the form

(1.12) ∞ n=1 g(n) 2n n 2 H n
in a convenient way in terms of a relatively "manageable" definite integral over an elementary function. As discussed in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], this is very useful because it is not obvious in general how to symbolically compute series of the form noted in (1.12) by applying parameter derivatives to known hypergeometric identities.

We remark that throughout the course of our present article, expressions such as "closed form" are meant to include evaluations involving the gamma function. Also, it is useful to note that it is convenient for our purposes to use the "Mathematica definition" of the complete elliptic integral of the first kind K whereby

K(k) = π 2 • 2 F 1 1 2 , 1 2 1 k ,
with E defined so that

E(k) = π 2 • 2 F 1 -1 2 , 1 2 1 k .

Motivating examples.

As we previously noted, the following theorem follows immediately from the formula for (1.11) given in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF]. We offer a new proof of this result, to illustrate how the main technique in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] can be used to evaluate series as in (1.8).

Theorem 2.1.

∞ n=1 2n n 2 Hn 32 n = Γ 2 ( 1 4 ) 4 √ π 1 -4 ln(2) π [13].
Proof. Letting f (n) be equal to 2 -n (2n -1) in Lemma 1.1, we find that the series given in the above theorem is equal to:

- 2 √ 2 π 1 0 ln 1 -x 2 √ 1 -x 2 √ 2 -x 2 dx - Γ 1 4 2 ln(2) π 3/2 .
Using the substitution u = 1 -x2 in the above integrand, we find that

∞ n=1 2n n 2 H n 32 n = - √ 2 π 1 0 ln(u) √ 1 -u 2 √ u du - Γ 1 4 2 ln(2) π 3/2 .
The Mathematica computer algebra system is able to evaluate the above integral directly, yielding the desired result. Alternatively, the substitution of the Maclaurin series for 11

To show how Lemma 1.1 may be applied in a non-trivial way to evaluate classes of variants of the summation given in Theorem 2.1, we consider the problem of evaluating the following natural analogue of the harmonic summation given in the above theorem.

(2.1)

∞ n=1 2n n 2 H n 32 n (n + 1)
It is not obvious as to how to evaluate this series following the techniques given in [START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF]. Mathematica 11 is not able to provide any kind of closed-form evaluation for the series in (2.1), and it is not obvious as to how to apply known integral formulas for harmonic numbers to evaluate this sum. For example, through an application of the formula

1 0 1 -x n 1 -x dx = H n
we see that the summation in (2.1) may be expressed as

1 0 8E( x 2 )+4(x-2)K( x 2 ) πx - √ 2 π Γ( 3 4
)

Γ( 5 4 ) x -1
dx, but it is not at all obvious how to compute the difficult integral given above. Similarly, through the use of the integral identity whereby

H n = -n • x∈[0,1] ln(1 -x) • x n-1 dx
for n ∈ N 0 , we find that the problem of symbolically computing (2.1) is equivalent to the difficult problem of evaluating the following integrals.

-1 16

1 0 2 F 1 3 2 , 3 2 3 x 2 ln(1 -x) dx = 4 π 1 0 E x 2 x 2 + x 2 -2 K x 2 x 2 ln(1 -x) dx
In a similar fashion, it appears that it would be infeasible to make use of known integral formulas for central binomial coefficients or Catalan-type numbers to determine a closed-form expression for the infinite series in (2.1). To illustrate this assertion, if we factor out the expression C n = 1 n+1 2n n in the summand in this series and replace this factor with a standard integral expression for the Catalan number C n of order n, we see that the infinite sum in (2.1) may be expressed as 2

√ 2 π 4 0 4 -x x(8 -x) ln 1 2 + 2 8 -x dx,
which cannot be evaluated by state-of-the-art CAS programs. Similarly, if we substitute a Wallis-type integral into the summand in (2.1), this would yield a very recalcitrant integral such as that given below.

1 π 2π 0 4 sec 2 (t) 1 - cos 2 (t) 2 ln 2 1 - cos 2 (t) 2 - 1 - cos 2 (t) 2 + 1 ln 1 - cos 2 (t) 2 + 1 + ln(2) dt
The above discussion shows that it is not feasible to use standard or conventional integration methods to evaluate

∞ n=1 2n n 2 H n 32 n (n + 1) ,
which shows, in part, why determining the closed-form evaluation given below is challenging. Our proof of the following theorem nicely illustrates how the technique introduced in [4] can be applied in a creative way to produce a simple closed-form evaluation that does not follow immediately from Lemma 1.1.

Theorem 2.2. The series

∞ n=1 1 32 n 2n n 2 H n n + 1 is equal to 8 - 2Γ 1 4 2 π 3/2 - 4π 3/2 + 16 √ π ln(2) Γ 1 4 2 . 13 Proof. Letting f (n) = 1 2 n 2n-1
n+1 , by Lemma 1.1, we have that the series

∞ n=0 1 32 n ( 2n n )
2 Hn n+1 is equal to the following expression:

4 π 1 0 √ 4 -2x 2 -2 ln 1 -x 2 x 2 √ 1 -x 2 dx - 16 √ π ln(2) Γ 1 4 2 .
So, the problem of computing the series given in Theorem 2.2 is equivalent to the problem of evaluating the following integrals, neither of which can be directly evaluated in closed form by state-of-the-art computer algebra systems.

1 0 √ 4 -2x 2 -2 ln 1 -x 2 x 2 √ 1 -x 2 dx = 1 0 √ 2 √ 1 + u -2 ln(u) 2(1 -u) 3/2 √ u du
Rewriting the expression

1 0 -2 + √ 2 √ 1 + u ln(u) (1 -u) 3/2 √ u du as (2.2) 4π + √ 2 1 0 √ 1 + u ln(u) (1 -u) 3/2 √ u du,
we find that current CAS software also cannot compute this latter integral.

Our strategy for computing (2.2) in closed form is to: find a formula for the Maclaurin series coefficients for the expression √ 1+u

(1-u) 3/2 in the integrand in (2.2), then multiply each term in the corresponding series expansion by ln(u) √ u , and then integrate term-by-term. This may appear to be a very roundabout way of determining the symbolic value of the definite integral in (2.2), but it is not at all clear what kinds of integration methods could be successfully applied to find (2.2) in closed form.

Since (2.3)

√ 1 + u (1 -u) 3/2 = ∞ n=0 u 2n (1 + 2u + 4n(1 + u))Γ 1 2 + n √ πΓ(1 + n) we have that √ 1 + u (1 -u) 3/2 = ∞ n=0 (1 + 4n)u 2n Γ 1 2 + n √ πΓ(1 + n) + ∞ n=0 4u 1+2n Γ 3 2 + n √ πΓ(1 + n)
so that we may evaluate the definite integral in (2.2) explicitly as follows:

1 0 √ 1 + u ln(u) (1 -u) 3/2 √ u du = ∞ n=0 (1 + 4n)(-4)Γ 1 2 + n (1 + 4n) 2 ( √ πΓ(1 + n)) + ∞ n=0 4(-4)Γ 3 2 + n (3 + 4n) 2 ( √ πΓ(1 + n)) = - Γ 1 4 2 + πΓ 3 4 2 √ 2π .
We thus obtain the desired result.

As discussed in the Introduction, the integration strategy given in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] is much more versatile compared to the use of specific hypergeometric identities to produce specific results on harmonic summations. To further illustrate this idea, we offer a complete proof of the following new evaluation. Proof. We observe that

∞ n=0 32 -n 2n n 2 H n 2n -1 is equal to (2.4) 2 1 0 √ 2-x 2 ln(1-x 2 ) √ 2 √ 1-x 2 dx + π ln(2) 2 √ π Γ( 1 4 ) 2 + 4Γ( 5 4 ) 2 π 3/2 π NEW SERIES INVOLVING 2n n 2 Hn
from the main lemma given in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], and we thus observe that the series given in Theorem 2.3 is also equal to (2.5)

2 1 0 √ 1+u ln(u) √ 1-u √ u du 2 √ 2 + π ln(2) 2 √ π Γ( 1 4 ) 2 + 4Γ( 5 4 ) 2 π 3/2 π .
The Mathematica computer algebra system is able to directly evaluate the datum from (2.5), yielding the desired result.

Theorem 2.4. The series

∞ n=1 1 32 n 2n n 2 H n (2n -1) 2
is equal to

(π -4 ln(2))Γ 1 4 2 8π 3/2 - 2 π (π + 4 ln(2) -6)Γ 3 4 Γ 1 4 .
Proof. We begin by noting that a direct application of the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] shows that the series given in the above theorem, which cannot be evaluated using computer algebra systems such as Mathematica 11, is equal to the following expression.

- √ 2 π 1 0 √ 2 -x 2 + x sin -1 x √ 2 ln 1 -x 2 √ 1 -x 2 dx - 2 π Γ 3 4 2 + 2Γ 5 4 2 ln(2) Γ 3 4 Γ 5 4
Expanding the above integrand, we again encounter the integral

1 0 √ 2 -x 2 ln 1 -x 2 √ 1 -x 2 dx
which we had previously seen in (2.4). So, we find that the infinite sum given in Theorem 2.4 is also equal to the following.

- √ 2 π 1 0 x sin -1 x √ + Γ 1 4 (π -4 ln(2)) 4 √ 2πΓ 3 4 - Γ 3 4 2 (π + 8 ln(2) -4) 2π 3/2
By rewriting the above integral as (2.6)

1 0 sin -1 √ 1-u √ 2 ln(u) 2 √ u du
and then substituting the expression

√ 1-u √ 2
into the Maclaurin series for the inverse sine, we see that the definite integral given in (2.6) is also equal to

-2 -5 2 π ∞ n=0 2 -5n 2n n 2 H n+1 n + 1 - √ πΓ 3
4 ln(4) 4Γ 5 4 and we may thus apply Theorem 2.2 to yield the desired result.

3. Generalizations and variations. We begin by noting that one of the key ingredients in our proof of Theorem 2.2 was based on the the use of a Maclaurin-type series for the expression (3.1)

√ 1 + u (1 -u) 3/2
for u ∈ [0, 1), and that the determination of a suitable power series expansion for this expression was nontrivial in that it is not obvious as to how to find explicit formulas for the Taylor series coefficients for (3.1) without already knowing the formula we had introduced in (2.3). For our strategy in computing the Maclaurin series coefficients for (3.1), we had planned to make use of known results on the series expansion for 1 + u 1 -u and then "interpret" the left-hand factor in the left-hand side of

1 1 -u 1 + u 1 -u = √ 1 + u (1 -u) 3/2
as a partial sum operator. The preceding discussion nicely illustrates how successful applications of the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] often require creative manipulations of generating functions.

It seems that applying the proof technique for Theorem 2.2 to try to determine the value of

(3.2) ∞ n=1 2n n 2 H n 32 n (n + 2)
may be cumbersome or infeasible. However, by means of a simple re-indexing argument, we see that there is a very elegant connection between (3.2) and Theorem 2.2, Theorem 2.3, and Theorem 2.4.

Theorem 3.1. The series ∞ n=1 2n n 2 H n 32 n (n + 2) is equal to 64 9 - Γ 1 4 2 (π -4 ln(2) + 18) 9π 3/2 - 2 2 π Γ 3 4 (9π + 36 ln(2) -16) 9Γ 1 4 .
Proof. Apply the re-indexing technique outlined below.

∞ n=1 2n n 2 H n 32 n (n + 2) = ∞ n=2 1 32 n-1 2n-2 n-1 2 H n-1 n + 1 = ∞ n=2 1 32 n-1 2n-2 n-1 2 H n -1 n n + 1 = ∞ n=2   32 1-n 2n-2 n-1 2 H n n + 1 - 32 1-n 2n-2 n-1 2 n(n + 1)   = - 2 π Γ 1 4 -3Γ 7 4 6Γ 7 4 + 32 ∞ n=2 32 -n 2n-2 n-1 2 H n n + 1 = - 2 π Γ 1 4 -3Γ 7 4 6Γ 7 4 + 8 ∞ n=2 2 -5n n 2 2n n 2 H n (n + 1)(2n -1) 2 = - 2 π Γ 1 4 -3Γ 7 4 6Γ 7 4 + 8 ∞ n=2 2 -5n 2n n 2 H n 1 9(n + 1) + 1 6(2n -1) 2 + 5 18(2n -1)
We now have that the desired result follows immediately from Theorem 2.2, Theorem 2.3, and Theorem 2.4.

To compute a given series of the form

∞ n=1 2n n 2 H n 32 n (n + z)
in closed form for z ∈ Z >0 , we make use of the inductive approach described below. We begin by re-writing the above summation as suggested below.

∞ n=1 2n n 2 H n 32 n (n + z) = ∞ n=1 2n n 2 H n 32 n ((n + 1) + z -1) = ∞ n=2 2n-2 n-1 2 H n-1 32 n-1 (n + z -1) = 8 ∞ n=2 2 -5n n 2 2n n 2 H n-1 (2n -1) 2 (n + z -1) = 8 ∞ n=2 2 -5n n 2 2n n 2 H n -1 n (2n -1) 2 (n + z -1) = 8 ∞ n=2 2 -5n n 2 2n n 2 H n (2n -1) 2 (n + z -1) -8 ∞ n=2 2 -5n n 2n n 2 (2n -1) 2 (n + z -1)
It is easily seen that hypergeometric series of the form

n 1 32 n n 2n n 2 (2n -1) 2 (n + m) NEW SERIES INVOLVING 2n n 2 Hn
19 always have closed-form evaluations for m ∈ N, as may be verified by writing

∞ n=0 n 2n n 2 x n+m-1 32 n (2n -1) 2 = 1 8 x m 2 F 1 1 2 , 1 2 2
x 2 and evaluating the above expressions as

1 8 x m 4 1 -2 x K x 2 π + 8E x 2 πx
and then using known results on moments of complete elliptic integrals. Now, we observe that we may expand the factor

n 2 (2n -1) 2 (n + z -1)
from the summand in the series

n n 2 (2n -1) 2 (n + z -1) 1 32 n 2n n 2 H n as follows. n 2 (2n -1) 2 (n + z -1) = z -1 2z -1 2 • 1 n + z -1 + 4z -3 2(2z -1) 2 • 1 2n -1 + 1 2(2z -1) • 1 (2n -1) 2
So, in our attempts to compute

∞ n=1 2n n 2 H n 32 n (n + z)
in closed form, we see that this problem amounts to the symbolic computation of the following series.

∞ n=1 2n n 2 H n 32 n (n + (z -1)) ∞ n=1 2n n 2 H n 32 n (2n -1) ∞ n=1 2n n 2 H n 32 n (2n -1) 2
So, we see that Theorem 2.2 may be regarded as the "base case" for our inductive technique, with Theorem 2.3 and Theorem 2.4 providing the required evaluations for the latter two sums given above, thus highlighting the import and the utility of the theorems given in Section 2. Through an application of the technique described above, we obtain the following new results. We encounter computational obstacles in attempting to apply the same kind of inductive procedure with respect to summations of the form n∈N 2n n 2 H n 32 n (2n -2z + 1) 21 for z ∈ Z >0 . This illustrates how the problem of computing generalizations of (1.7) can be very complicated and often requires a degree of ingenuity in the application of the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF]. The problem of evaluating

∞ n=1 2n n 2 H n 32 n (n + 3) = 2048 225 - 4 √ π(297π + 1188 ln(2) -800) 225Γ 1 4 2 - 2Γ 1 4 2 (313 + 25π -100 ln(2)) 225π 3/2 ∞ n=1 2n n 2 H n 32 n (n + 4) = 16384 1225 + Γ 1 4 2 (-15974 -1425π + 5700 ln(2)) 3675π 3/2 - 4 √ π(-7984 + 2401π + 9604 ln (2)) 1225Γ 1 4 2
(3.3) ∞ n=0 2n n 2 H n 32 n (2n -3)
is especially interesting, since there is a beautifully unexpected corollary of the symbolic evaluation for the infinite summation given above that we discuss in Section 5.

Suppose that we were to attempt to evaluate the series in (3.3) using Theorem 2.3, following the inductive strategy that had been employed in our generalization of the proof for Theorem 3.1. So, we must re-index (3.3) as demonstrated below. (n + 1) 3 (2n -1) in closed form. Moreover, if we apply partial fraction decomposition with respect to the factor (2n + 1) 2 (n + 1) 2 (2n -1) in the "re-indexed" harmonic summation given above, we obtain the expression 16 9(2n -1) + 10 9(n + 1) -1 3(n + 1) 2 , but it is not at all clear as to how the series n∈N 2n n 2 H n 32 n (n + 1) 2 could be evaluated, even through an application of the main lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], since this would require the symbolic computation of a difficult integral such as

∞ n=1 32 -n 2n n 2 H n 2n -3 ∞ n=0 32 -n-1 2n+2 n+1 2 H n+1 2n -1 1 8 ∞ n=0 2 -5n (2n + 1) 2 2n n 2 H n + 1 n+1 (n + 1) 2 (2n -1) 1 8 ∞ n=0 2 -5n (2n + 1) 2 2n
1 0 ln(u) ln 1 + √ 1+u √ 2
(1 -u) 3/2 √ u and it is not clear as to how to apply the Maclaurin series substitution strategy employed to evaluate (2.2).

To evaluate series of the form

n∈N 2n n 2 H n 32 n (2n -2z + 1)
for z ∈ Z >0 , we may apply the following procedure inspired by our proofs for Theorem 2.2 and Theorem 2.3. We leave it as an exercise to verify the efficacy of this algorithm.

(1) Let f (n) = 2n-1 2 n (2n-9) , apply the funamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], and evaluate the corresponding hypergeometric series (2) Through the application of Lemma 1.1 noted above, we obtain an integrand of the form

2 -x 2 1 -x 2 ln 1 -x 2 p(x)
for a polynomial p(x) with algebraic coefficients. Apply the substitution u = 1 -x 2 ; in closed form for z ∈ Z >0 , and we leave it as an exercise to formalize this idea. For example, using an analogue of the procedure described above, we obtain the following result. for z ∈ Z >0 had not been discussed in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF]; series of the form in (4.1) cannot be evaluated through a direct or straightforward application of Lemma 1.1, as is easily verified. We thus make use of a more "unusual" approach to prove the following result, which serves as the base case for an inductive generalization. Proof. We begin by making use of the following result that had been introduced in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF]:

∞ n=1 2n n 2 H n 16 n (2n -3) = -68 + 120 ln(2) 27π .
We remark that the above result follows directly from the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] by letting f (n) = 2n-1 2n-3 . Now apply the following reindexing argument.

-68 We remark that for the final equality given above, we implicitly made use of Lemma 1.1 in the case whereby f (n) = 1. So, we have shown that the problem of symbolically computing the infinite series given in the above theorem is equivalent to the remarkably simpler problem of computing the hypergeometric series given below. Surprisingly, Mathematica is not able to compute the above series. Making use of the Catalan number integral formula whereby

+ 120 ln(2) 27π = ∞ n=1 16 -n 2n n 2 H n 2n -3 = ∞ n=0 16 -(n+1) 2n+2 n+1 2 H n+1 2n -1 = 1 4 ∞ n=0 4 -2n (2n + 1) 2 2n n 2 H n+1 (n + 1) 2 (2n -1) = 1 4 ∞ n=0 4 -2n 2n n 2 H n - 1
2n n n + 1 = 1 2π 4 0 x n 4 -x x dx,
we find that the series in (4.2) is also equal to

1 2π 4 0 -8 √ 4 -x 1 x 3/2 -2 + √ 4 -x+ 2 ln(2) -2 ln 1 + 1 - x 4 dx.
which Mathematica is able to evaluate as -32G π -16 + 48 π + 16 ln(2), thus completing our proof.

By analogy with our generalization of Theorem 3.1, we make use of the following procedure to compute infinite series of the form

∞ n=0 2n n 2 H n 16 n (n + z) 2
for z ∈ Z >0 . Begin by rewriting (4.1) as follows.

∞ n=1 2n n 2 H n 16 n (n + z) 2 = ∞ n=2 2n-2 n-1 2 H n-1 16 n-1 (n + z -1) 2 = 4 ∞ n=2 4 -2n n 2 2n n 2 H n (2n -1) 2 (n + z -1) 2 -4 ∞ n=2 4 -2n n 2n n 2 (2n -1) 2 (n + z -1) 2 It is easily seen that n 4 -2n n 2n n 2 (2n -1) 2 (n + z -1) 2
has a closed-form expression for z ∈ Z >0 , and we leave it as an exercise to show this. Through the use of partial fraction decomposition, expand the quotient n 2 (2n -1) 2 (n + z -1) 2 as suggested below.

-2z -2 (2z -1) 3 •

1 n + z -1 + 1 (2z -1)
This is very interesting because it is unexpected that the main integration technique from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] can also be applied to determine new series involving alternating harmonic numbers, and this suggests that there may be a deep connection between the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] and the main integration technique from [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF]. How can we obtain new classes of series containing harmonic-like numbers of the form H 2n as in (5.2) by showing how the definite integral in Lemma 1.1 can be expressed in terms of parameter derivatives of hypergeometric expressions of the form p F q (-1), as above?

16 n

 16 (2n -1) 2

  the derivation of one of the main identities given in[START_REF] Guillera | More hypergeometric identities related to Ramanujan-type series[END_REF].

Theorem 2 . 3 .

 23 The following equality holds:

n∈N 2n n 2 32

 2 n (2n -2z + 1) in closed form;

Theorem 4 . 1 .

 41 The following equation holds:

  4. Ramanujan-type formulas. Inspired by the method we had applied to generalize Theorem 3.1, we consider the use of similar strategies to determine explicit evaluations for new1 32 n which served as something of a basis for Section 2 and Section 3. To illustrate this idea, we begin by considering the problem of finding a symbolic evaluation for the simple and natural-looking series

				π series with
	summands containing factors of the form
			2n n 16 n , 2 H n
	which were the main kinds of mathematical objects under investigation
	in [4], as opposed to expressions of the form
			2n n	2 H n
		∞ n=1	C 2 n H n 16 n ,
	letting C n denote the n th Catalan number, as above. Summations of
	the form
	(4.1)	∞ n=1	2n n 16 n (n + z) 2 2 H n

  Starting with Theorem 4.1 as the base case, we may thus evaluate (4.1) recursively since it is easily seen that series of the form16 n (n + m)for m ∈ N may be evaluated directly through Lemma 1.1. Following this algorithm, we obtain the new results indicated below.as a factor lead us to a very surprising discovery concerning a new series involving alternating harmonic numbers, as elaborated below. It seems that the problem of generalizing this discovery would form a promising area of research.By Lemma 1.1, by letting f (n) = 2 -n (2n-1)Mathematica 11 is not able to evaluate the above expression in closed form. However, interestingly, if the function Simplify is applied to this expression, Mathematica produces the following output, after 11 minutes of computation, letting γ denote the Euler-Mascheroni constant, and letting regularized hypergeometric functions be denoted It can be shown that from the above output, together with the symbolic calculation for (5.1) that we had provided in Section 3, we have that the series(5.2)

				NEW SERIES INVOLVING 2n n	2 Hn	29
	= with p Fq . 6272 3375 ∞ n=1 2n -2 H n 1024 ln(2) 225 1 18 √ 2π 10 √ πΓ 1 4 ln(2) 23632 + 30720G -54208 ln(2) + 3375π Γ 3 4 + 24 √ πΓ 3 4 ln(2) + Γ 1 4 n 16 n (n + 4) 2 = 129536 128625 -4096 ln(2) 1225 783408 + 860160G -1484864 ln(2) 3 4γE(-1) + π -√ π 2 ∂ ∂x 3 F2 -1 2 , 1 2 , 1 2 -1 + 1, x x= 1 2 + 128625π 5. Conclusion. Our explorations on the symbolic computation of ∂ ∂x 3 F2 -1 2 , 1 2 , 1 2 2 x= 1 2, x -1 +
	2 • series with summands containing 1 32 2 ∂ ∂x 2 F 1 -1 2 , 1 2 -1 + x x=1 n 2n (2n -1) 2 1 2 n H n ∂ ∂x 2 F1 -1 2 , 1 2 -1 x x=2
						+	(z -1) 2 (2z -1) 2 •	1 (n + z -1) 2
	series			∞ n=1	+ -	4z -4 (2z -1) 3 • 1 n 2n 16 n	1 2n -1 2n-3 2 (2n + 3)H 2n , we find that the infinite (2n -1)(n + 1)
	∞ n=1 ∞ 32 -n 2n 2n n 2 H n n 2 H n n=0 2n -3 2 (π -4 ln(2)) is may be evaluated as (5.1) is equal to the following. 5Γ 1 4 12 √ 2π 3/2 -Γ 3 4 (3π + 12 ln(2) -20) 3 √ 4 πΓ 1	-
		2 π	0	1	√	2 -x 2 + x 2 √ 3 √ 2 √ 2 -x 2 ln 1 -x 2 1 -x 2	dx+
	∞ n=1 = 112 16 n (n + 2) 2 2n n 2 H n 27 -64 ln(2) 9 5Γ 1 4 ln(2) + 9 √ 2πΓ 3 4	27π 16(13 + 24G -44 ln(2)) 2 2 π Γ 3 4 ln(2) + 3Γ 1 4
	∞	2n n	2 H n
	n=1	16 n (n + 3) 2

√ 1-u

in the above integrand also may be used to symbolically compute the above integral.

ln 1 -x 2 √ 1 -x 2 dx

23

(3) We thus obtain an integrand of the form

for a polynomial q(u) with algebraic coefficients. Replace the expression 1+u 1-u with its Maclaurin series

in the above integrand and integrate term-by-term.

Using the above procedure, we obtain the following results, thus illustrating the versatility of the main lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF].


Alvaro Javier Fuentes Suárez
email: alvaro.fuentes-suarez@inria.fr

Evelyne Hubert

Scaffolding skeletons using spherical Voronoi diagrams: feasibility, regularity and symmetry

Keywords: feasibility, regularity and symmetry Quad mesh generation, Skeleton based modeling, Procedural modeling, Scaffold

published or not. The documents may come

Introduction

Skeletons are used in 3D graphics for modeling and animating articulated shapes. The user can design a complex shape by sketching a simple geometric object that is the input to a surface generating algorithm. By making changes in the skeleton it is possible to change the shape in an intuitive way. Due to their low dimensional nature, skeletons can serve as an efficient and compact representation of a surface. In this context a skeleton-based mesh generation method is needed.

The idea developed here is to construct a "coarse" quad mesh that tightly follows the structure of the skeleton. Following the terminology of [START_REF] Panotopoulou | Scaffolding a skeleton[END_REF] we call this process scaffolding and the corresponding coarse mesh a scaffold. This is used as an intermediate step in many applications. A scaffold can be used to generate a surface, either by subdivision as in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF], or as initial patchwork for a spline surface [START_REF] Blidia | G1-smooth splines on quad meshes with 4-split macro-patch elements[END_REF][START_REF] Karčiauskas | Curvature continuous bi-4 constructions for scaffold-and sphere-like surfaces[END_REF] (see Figure 1). It is used as an intermediate step in the extraction of a quad layout on a given triangular mesh [START_REF] Baerentzen | Interactive shape modeling using a skeleton-mesh co-representation[END_REF][START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF], and for compatible quadrangulation [START_REF] Yao | Compatible quadrangulation by sketching[END_REF]. An application we present here is the polygonization of skeleton-based implicit surfaces into quad-dominant meshes. In particular, we use our method for visualization of convolution surfaces [START_REF] Bloomenthal | Convolution surfaces[END_REF][START_REF] Zanni | Skeleton-based Implicit Modeling & Applications[END_REF]. Our main contribution is in the theoretical foundations of an algorithm that computes a scaffold for any skeleton, independently of its topology, and with no user interaction.

In this paper we deal with skeletons made of line segments that do not intersect except at the endpoints, then called joints.

Previous work

One of the earliest ideas for a scaffold construction was to sweep a fixed polygonal cross profile along the segments, stitching the generated quads by means of a convex hull construction at the joints. This idea was first used in [START_REF] Srinivasan | Solidifying wireframes[END_REF], producing meshes with some triangular faces resulting from the stitching process. For quadrilateral Figure 1: Recreation of some of the scaffolds used in [START_REF] Karčiauskas | Curvature continuous bi-4 constructions for scaffold-and sphere-like surfaces[END_REF], automatically computed by our method from suitable skeletons.

cross profiles B-meshes method [START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF] improved upon this by merging triangles in order to get a quad-dominant mesh. Still the stitching process might leave some triangular patches at the joints. [START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF] and [START_REF] Yao | Compatible quadrangulation by sketching[END_REF] proposed a scaffolding technique that generates a pure-quad mesh by extruding boxes emanating from cubes positioned at the joints. In [START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF] the subdivision of the cubes is modeled with an Integer Linear Program, for which a solution might fail to exist when there are cycles in the skeleton. "Lids" [START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF]Figure 6] are introduced as a workaround. [START_REF] Yao | Compatible quadrangulation by sketching[END_REF] also introduces extra quads around joints [START_REF] Yao | Compatible quadrangulation by sketching[END_REF]Figure 4].

The method we propose is based on the Skeleton to Quaddominant polygonal Mesh (SQM) method in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF] which was limited to skeletons without cycles. SQM first defines the mesh points around the joints and then recreate a quad based "tubular" polyhedral surface around each line segment. SQM can be regarded as a three steps process:

1. Partition the unit sphere centered at the joints into regions, one for each incident line segment.

2. Discretize each region into a cell (as an ordered set of points on its boundary) such that the two cells at the extremities of each line segment are compatible, i.e. have equal number of points. Points on the boundary of two regions are part of the two corresponding cells.

3. Link, in a bijective way, the cells at the extremities of each line segment. These links define the quads on the mesh.

For Step 3, SQM [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF] defines the links by minimizing the total length of the line segments they define. In Step 2 Baerentzen et al. propose an algorithm for inserting additional vertices on the cells in such a way that the compatibility constraint is satisfied. Yet this algorithm does not work in the presence of cycles [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]. The existence of a possible discretization was actually not proved. Furthermore there was no analysis on the optimality of SQM with respect to the number of quads in the scaffold.

The partition of the sphere, in Step 1, used in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF] can be recognized to be a Voronoi diagram on the sphere. [START_REF] Panotopoulou | Scaffolding a skeleton[END_REF] introduces a partition of the sphere in quadrangles that makes the compatibility of cells (Step 2) trivial. Yet the partition is not canonical and the convexity of the regions is not guaranteed.

Contributions

Our method follows SQM [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF] but we address and solve the crux difficulty that represents Step 2 for skeletons of arbitrary topology. We formalize the creation of compatible cells (Step 2) as an Integer Linear Program (IP) that minimizes the total number of quads. We prove that, for a Voronoi partition of the spheres, there exists a solution for the IP (feasibility), even in the presence of cycles. There is thus always an optimal solution that can be computed by an IP solver. Feasibility is proven thanks to a numerical characterization by Rivin [START_REF] Rivin | A characterization of ideal polyhedra in hyperbolic 3-space[END_REF] of graphs combinatorially equivalent to inscribable polyhedrons (i.e. with vertices on a sphere) that applies to the dual of the Voronoi diagram. Our minimization criterion is what ensures the coarsest mesh among those based on Voronoi partition of the spheres and additional pragmatic geometric restrictions. The solution of the IP determines the cross profile on each segment.

We present two other constructions to generate symmetric and regular scaffolds. In the former case the scaffold respects the symmetries of the skeleton. In the latter case, the scaffold has equal number of quads around each line segment of the skeleton, ensuring a similar cross profile for each line segment. Both possibilities are natural requirements for geometric modeling. With either or both requirements, we prove the feasibility of the constraints and are thus in a position to compute optimal solutions in the total number of quads.

To the extent of the knowledge of the authors the only paper that integrates the symmetries of the skeleton into the computation of the scaffold is [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF], with the limitation that only one reflection symmetry can be taken into account for each skeleton (in addition to being restricted to cycle-free skeletons). Here we present a much more general approach that is able to compute scaffolds that respect any group of symmetries of the skeleton.

The article is organized as follows. In Section 2 we formalize the notions of skeleton and scaffold. In Section 3 we prove the existence of: standard, regular, and symmetric scaffolds, for any topology. In Section 4 we define an objective function and introduce the IPs that find optimal solutions. The algorithms to construct a scaffold are detailed in Section 5 with some further discussion in Section 6. An application to polygonization of skeleton-based implicit surfaces is shortly presented in Section 7.

A sketch of the first feasibility proof was presented at the conference LAGOS 2017. An extended abstract of the talk given there is available in the proceedings [START_REF] Fuentes Suárez | Scaffolding skeletons using spherical voronoi diagrams[END_REF]. Here we generalize and extend the result in [START_REF] Fuentes Suárez | Scaffolding skeletons using spherical voronoi diagrams[END_REF]. The precise objective function of the IP, the details of the algorithms, as well as the regular and symmetric cases are presented in this paper for the very first time.

Skeletons & Scaffolds

In this paper a skeleton is a finite set S of spatial line segments satisfying the following property: any two line segments intersect at most at one of their endpoints. A skeleton S defines naturally a graph G S = (V S ; E S ) by identifying the set of nodes V S with the set of all endpoints in S, and the set of edges E S with the set S itself. An edge e P E S connecting two nodes a; b P V S can be alternatively represented as ab. If e = ab we say that e is incident to a (or b) and write e (a (or e (b). A node with only one incident edge is called dangling node while a node with more than one incident edge is called a joint. A node with precisely two incident edges is called an articulation. We denote by L S & V S the set of all dangling nodes, M S & V S the set of all articulations, and N S = V S (L S M S ) the set of the remaining joints.

The sphere centered at v P V S with radius ε v > 0 is denoted S v , and A v = feS v j e P E S ; e (vg is the set of the points that are the intersection of the line segments incident to v with S v (Figure 2b). For what follows the choices of ε v (v P V S ) is independent of our method. We assume though that no two spheres intersect. Note that, once the scaffold mesh is created, max vPV S ε v is an upper bound to the distance between any point in the edges of the scaffold and the skeleton. end-points of the arcs delimiting R v e and some additional points (possibly none) per arc (Figure 2d). The points chosen in one arc define a polyline that represents the arc (Figure 2e). The number of segments in the polyline is called number of subdivisions of the arc, it is one less than the number of points taken in the arc.

One should also consider some geometrical constraints. The number of points in a cell must be at least 3 (4 is more customary [START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF][START_REF] Panotopoulou | Scaffolding a skeleton[END_REF][START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF][START_REF] Yao | Compatible quadrangulation by sketching[END_REF]). Long arcs, i.e. arcs with length greater than or close to π, must be subdivided into at least two segments. Examples of degenerate cases for cells with at least 3 or 4 points are shown in Figure 3. If the node v is an articulation, the arc separating the two Voronoi regions in A v is a circle. Thus an associated cell C v e consists of points on the boundary circle. For v a dangling node, Vor(A v ) consists of a single region with no boundary. In this case the associated cell C v e consists of points on the circle defined by the sphere S v and the plane through v normal to e. In both cases the points on the cells are given by only one arc (great circle) and are taken such that the planar polygon they define encloses v.

There seems to be a preference in the literature for at least quadrangular cross profiles [START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF][START_REF] Panotopoulou | Scaffolding a skeleton[END_REF][START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF][START_REF] Yao | Compatible quadrangulation by sketching[END_REF]. To guarantee this the cells must have at least four points. We show in Figure 14 that three is also an adequate choice. In general, arcs can be discretized as a single segment (i.e. with no additional point) but the extra restrictions on the minimum number of points on each cell must be enforced to avoid cases like in Figure 3.

A scaffold K S is defined as a pair (P S ; Φ S ), satisfying 1. P S = fC v j v P V S g, where each C v = fC v e j e P E S ; e ( vg is a family of cells representing a partition of S v according to Vor(A v ).

2. Φ S = fφ e j e P E S g is a family of bijections φ e between C a e and C b e for e = ab.

For e = ab we say that C a e and C b e are linked cells. Similarly, if C a e = hp 1 ; p 2 ; : : : ; p n i we say that p i is linked with φ e (p i ) and the pair hp i ; φ e (p i )i is called a link. The realization of a scaffold as a mesh is through the quads defined by the four-points tuples hp i ; φ e (p i ); φ e (p i+1 ); p i+1 i.

To construct the links we follow the same strategy as in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]: provided the cells have the same number of points, choose the bijection where the total length of the segments defined by the links is minimal.

The existence of the scaffold, for a given skeleton, is thus established if and only if we can discretize the Voronoi regions into compatible cells: we need the cells C a e and C b e (for all e = ab P E S ) to have the same number of points. This is far from obvious, specially in the presence of cycles in the skeleton.

Existence of scaffolds

We construct the set of linear equations over the integers that preside over the existence of a scaffold. The achievement is to prove the existence of a positive solution to this system, hence proving the existence of a scaffold for any given skeleton. The proof strongly relies on a property of Voronoi diagrams on the sphere, and more precisely on their duals. We examine this property first. We then prove the existence of a scaffold, which we qualify as standard. In geometric modeling it is desirable to have scaffolds that respect the symmetries of the underlying skeleton (symmetric scaffold), or to require a regularity on the number of quads around the line segments (regular scaffold). We can even seek scaffolds that satisfy both properties. We prove the existence of all these scaffolds.

Locally uniform discretization

The Delaunay triangulation Del(A v ) is the dual of Vor(A v ) [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF]. For each v P V S let E v be the set of edges of Del(A v ). Each edge in E v represents a common boundary between two regions in Vor(A v ).

For f P E v we define a positive integer x v f representing the number of subdivisions to be done to the corresponding arc (i.e. the number of segments in the polyline representation of the arc). For dangling nodes E v = / 0, we nonetheless introduce a phantom edge ν v , with an associated variable x v ν v , so that actually E v = fν v g.

The following lemma asserts that the Voronoi regions can be discretized uniformly, i.e. with an equal number of points. It is our main ingredient in proving the existence of scaffolds. Lemma 3.1. For v P V S , the local linear system

∑ f PE v f ((S v e) x v f = λ v Ve (v;ePE S (1)
has a solution ( xv f ; λv ) with positive integer entries.

We denote the solution ( xv f ; λv ) as local solution associated to the local system (1).

To guarantee a positive solution we rely in a proposition due to Rivin [START_REF] Rivin | A characterization of ideal polyhedra in hyperbolic 3-space[END_REF] (statement extracted from [START_REF] Dillencourt | Graph-theoretical conditions for inscribability and delaunay realizability[END_REF]) that gives a numerical characterization for a graph of inscribable type (i.e. a graph combinatorially equivalent to a polyhedron inscribed on a sphere [START_REF] Dillencourt | Graph-theoretical conditions for inscribability and delaunay realizability[END_REF]). Proposition 3.2 (I. Rivin). If a graph is of inscribable type then weights w can be assigned to its edges such that:

(i) For each edge e, 0 < w(e) < 1=2.

(ii) For each vertex v, the total weight of all edges incident to v is equal to 1.

Proposition 3.2 applies to Del(A v ) that is combinatorially equivalent to the convex hull of A v [START_REF] Brown | Geometric Transforms for Fast Geometric Algorithms[END_REF][START_REF] Grima | Computational Geometry on Surfaces[END_REF], and hence of inscribable type. It thus guarantees a positive real solution for [START_REF] Angelidis | Adaptive implicit modeling using subdivision curves and surfaces as skeletons[END_REF]. Note that guaranteeing a positive solution for a linear system is not a trivial task.

The following claim then relates the existence of an integer positive solution of a homogeneous linear system to the existence of a real positive solution.

Proposition 3.3. A homogeneous linear system with integer coefficients has a positive integer solution whenever it has a positive real solution.

of Proposition 3.3. Let Ay = 0 be a homogeneous linear system where A is a n ¢m matrix with integer entries. Observe that if there is a rational solution p with p P Q m , we can get an integer positive solution multiplying p by the least common multiple of denominators in the entries of p. Thus it is enough to prove that the system has a rational positive solution.

Let ỹ P R m be the real solution of the homogeneous system with positive entries, this implies that the set of solutions of the system is a (non-trivial) subspace. Since A has integer entries we get a rational basis for the solution space. Let fy 1 ; y 2 ; : : : ; y k g & Q m (k ! 1) be a basis of the solution space of A. We have that ỹ = ∑ k i=1 ci y i for some real coefficients ci . Let f : R k 3 R m be a function mapping

c P R k to f (c) = ∑ k i=1 y i c i . Let U = (0; ∞) m . Clearly U is open, and f is continuous, thus V = f 1 (U ) & R k is open. Moreover c = ( c1
; : : : ; ck ) P V hence V is not empty. The set of rational points in R k is dense, thus there exists q P Q k V. On the other hand f (q) P U, that is all the entries of f (q) are positive and rational and by definition f (q) is in the solution space of A. Therefore f (q) is a rational solution with positives entries for the system.

of Lemma 3.1. For v a dangling node or articulation, a solution is trivially found. If v is a joint with at least three incident edges, the existence of a real positive solution for the local system (1) comes from the fact that Del(A v ) is combinatorially equivalent to the convex hull of A v [START_REF] Brown | Geometric Transforms for Fast Geometric Algorithms[END_REF][START_REF] Grima | Computational Geometry on Surfaces[END_REF] which is an inscribed polyhedron. Thus Proposition 3.2 guarantees the existence of a real positive solution given by x v f = w( f ) and λ v = 1, the result follows from Proposition 3.3.

Standard scaffold

The number of points in a cell C v e (e P E S and e (v) is given by

jC v e j = ∑ f PE v f ((S v e) x v f : (2)
For each edge e = ab P E S , there is a bijection φ e P Φ S in the scaffold K S between the cells C a e and C b e . This is possible only if both cells have the same number of points, which gives the following compatibility equations ∑ hPE a h((S a e) (3) By definition x v f is a positive integer for all v P V S ; f P E v . In Section 2 we discussed some additional geometric constraints. They can be written in the form

( x v f P Z; x v f ! 1 Vv P V S ; f P E v Λ i (x v f ) ! s i i = 1; 2; ::: (4)
where Λ i (x v f ) are linear forms on the variables x v f with non-negative integer coefficients (not all zeros), and s i > 0 are integer constants. These can capture the requirements of having at least 3 (or 4) points on each cell, as well as subdividing long arcs into at least two segments. The existence proof works for any set of (Λ i ; s i ) with nonnegative coefficients. A practical realization of (4) for the geometric constraints commented in Section 2 and that should also serve as reference for the reader, is given by 8 > < > :

x v f P Z; x v f ! m(x v f ) Vv P V S ; f P E v ∑ f PE v f ((S v e) x v f ! c
Vv P V S ; e P E S ; e (v (5) where c is 4 or 3; m(x v f ) is 1 if the length of the arc associated to x v f is less than π δ , and 2 otherwise. The choice of a small constant δ P (0; π) determines the long arcs. Yet other constraints that may arise in applications, like subdividing specific arcs into a greater number of pieces, or requiring specific cells to have a greater number of points, can also be modeled in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]. We say that the global system defined by ( 3) is feasible if it has a solution satisfying [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]. Such a solution gives a way to discretize each region in the partition of the spheres into compatible cells, and thus allows to construct a scaffold. Proving the existence of a positive (integer or real) solution for a linear system is not a trivial task. The main result in our paper is the formal proof of feasibility of (3) which is stated in the following theorem.

Theorem 3.4. For any skeleton S, the linear system given by (3) has a solution with the entries x v f (v P V S ; f P E v ) satisfying the constraints given in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]. Therefore there exists a scaffold for S.

Proof. Using Lemma 3.1, for every v P V S , we take ( xv f , λv ) as a local solution satisfying (1). Then we take xv f = s λ λv xv f for all v P V S ; f P E v , where λ = ∏ uPV S λu and s = max i s i is an integer constant. Using xv f as subdivisions for the arcs we get that all the cells have the same number of points: s λ , it follows then that the equalities in (3) are trivially satisfied. The factors s λ = λv guarantee that the constraints in (4) are also satisfied. Thus xv f is a solution to (3) satisfying (4).

The discretization constructed in the proof above is far from optimal. Optimality is dealt with in Section 4 with the help of Integer Linear Programming.

Symmetric scaffold

It is a desirable property for a scaffold to respect the symmetries of the underlying skeleton [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]. As illustrated in Figure 4 and6, a standard scaffold need not satisfy this property. In this section we first define what is a valid symmetry for the skeleton. We then give the additional restrictions needed in order to obtain a scaffold that respects the symmetries of the skeleton. A skeleton symmetry of S is an isometry T : R 3 3 R 3 such that (i) T (v) P V S Vv P V S , and (ii) T (e) P E S Ve P E S .

If the radii ε v are different then in order to have a symmetric scaffold they must satisfy ε v = ε T (v) for all v P V S ; T P T S . We assume this is guaranteed for symmetric skeletons.

Conditions ((i)) and ((ii)) say that T keeps the set of nodes V S and edges E S (hence G S ) invariant, thus the whole skeleton is kept fixed: T (S) = S. Since T is an isometry and e = ab P E S is the line segment between the nodes a and b (including both), then T (e) is the edge (line segment) connecting T (a) and T (b).

With conditions ((i)) and ((ii)) we avoid cases, as illustrated in Figure 5, where there is a geometric symmetry for S that does not map elements of G S into elements of G S . Thus not all symmetries are skeleton symmetries.

A B O C

Figure 5: A skeleton with a false symmetry: the central symmetry with respect to O, maps the skeleton to itself when considered as a curve but not as a graph (B is not mapped to another node).

The set of skeleton symmetries T S forms a finite group under composition. For simplicity we denote T R as T R.

We say then that a scaffold K S = (P S ; Φ S ) respects the skeleton symmetry T P T S if

C T (v)
T (e) = T (C v e ) Vv P V S ; e P E S ; e (v;

(6) and φ T (e) = T φ e T 1 Ve P E S : [START_REF] Blidia | G1-smooth splines on quad meshes with 4-split macro-patch elements[END_REF] Since Voronoi diagrams depend only on the distance between points, it follows that

Vor(A T (v) ) = T (Vor(A v )): (8)
Notice that it is also true that A T

(v) = T (A v ) and E T (v) = T (E v ).
Hence Equation ( 6) can be achieved as soon as the subdivisions are done in a symmetric way for symmetric arcs. This is possible if

x v f = x T (v) T ( f ) VT P T S ; v P V S ; f P E v : (9)
To guarantee [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF] given [START_REF] Botsch | Polygon Mesh Processing[END_REF], it is sufficient to subdivide the arcs into equal length subdivisions since arc-length is preserved under isometries: the arc corresponding to f P E v (v P V S ) with arc-length θ is subdivided into x v f sub-arcs of length θ =x v f . Equation ( 7) means that the links defined by the bijections in Φ S define symmetric line segments, and hence symmetric quads.

Solutions of (3) satisfying ( 4) and the extra constraints given by ( 9), yield compatible cells that respect the skeleton symmetries. The existence of such a solution is denoted as feasibility of the symmetric scaffold, and it is established in Theorem 3.6.

The proof of Theorem 3.6 relies on the following lemma.

Lemma 3.5. Let xv f be a solution of the linear system given by (3) satisfying the constraints in (4). Then xv

f = ∑ T PT S xT(v) T ( f )
is also a solution of (3) satisfying (9) and (4).

Proof. Let T P T S . Since T S is a group we have T S = fRT j R P T S g. Hence xT(v)

T ( f ) = ∑ RPT S xRT(v) RT ( f ) = ∑ RPT S xR(v) R( f ) = xv f : (10)

Thus xv

f satisfies (9). We prove now that xv f is a solution of the linear system given by [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF]. For e = ab P E S we have T (E a ) = E T (a) because of the symmetric property of Voronoi diagrams [START_REF] Bloomenthal | Convolution surfaces[END_REF], and hence of its dual. Therefore

From ( 11), ( 12) and ( 13) we get

T (g) Ve P E S : (14)
Summing the latter equation over T P T S proves that xv f is a solution of (3). Theorem 3.6. For any skeleton S with a group of symmetries T S , the linear system given by (3) has a solution with the entries x v f (v P V S ; f P E v ) satisfying the constraints given in (4) and (9). Therefore there exists a symmetric scaffold for S. An example of a skeleton with a rotation symmetry is shown in Figure 6. If one is not interested in a solution that respects all the symmetries of the skeleton, one can restrict the set T S to be the group generated by a subset of the skeleton symmetries.

Regular symmetric scaffold

We call a scaffold regular if it has the same number of quads around each line segment. This is not an automatic property of standard scaffolds, as can be seen in figure 6b and7. To get regular scaffolds we need compatible cells that have the same number of points, this means that (3) must be replaced by

∑ f PE v f ((S v e)
x v f = λ Vv P V S ; e (v;e P E S [START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF] with λ an additional free variable. Notice that λ is independent of v P V S and so (15) implies (3). Solutions to the linear system [START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF] satisfying (4) are called regular solutions. Regularity ensures that all segments have a similar cross profile: a polygon with λ sides. If there is a regular solution we say that (15) is feasible, which implies the existence of a regular scaffold.

It is possible to get a scaffold that is at the same time regular and symmetric. We just need to ensure, besides (4), the extra constraints in [START_REF] Botsch | Polygon Mesh Processing[END_REF]. A regular scaffold need not be symmetric (Figure 8a). Conversely, a symmetric scaffold is not necessarily regular (Figure 8b).

The existence of regular symmetric scaffolds is established in the following theorem. It is sufficient to prove the feasibility of regular symmetric scaffolds to automatically get the feasibility of regular scaffolds. Indeed a regular scaffold can be regarded as a regular symmetric scaffold with T S the trivial group consisting of only the identity symmetry.

T ( f ) = = ∑ T PT S ∑ gPE T (v)
g((S T (v) T(e))

xT(v)

g = ∑ T PT S j ĈT(v)
T (e) j = s λ jT S j:

Symmetry is an important property for the scaffolds as can be appreciated in Figure 9 computed for one of the examples in [START_REF] Karčiauskas | Curvature continuous bi-4 constructions for scaffold-and sphere-like surfaces[END_REF]. The "correct" symmetric scaffold respecting all the symmetries of the skeleton (same as Figure 1a).

Optimal scaffolds

We have established that for any skeleton S admitting a group of symmetries T S , there exist a standard scaffold, a symmetric scaffold, a regular scaffold, and a regular symmetric scaffold. In this section we address the computation of optimal scaffolds. We model each scaffolding problem as an Integer Linear Program (IP) [START_REF] Chen | Applied Integer Programming[END_REF] that looks for a solution with a minimal number of quads. For this we need first to define an objective function, then the constraints of the IPs. The existence of a minimal solution is guaranteed by the existence theorems on Section 3, thus such minimal solutions can be computed with an IP solver.

Objective function

We want to express the total number of quads in a scaffold in terms of the subdivision of the arcs. For this we proceed as follows. Let Q be the set of all quads of a scaffold, and P the number of pairs hp;Qi such that p P Q P Q. Then we have that jQj = 1 4 P:

Indeed, the last equation is a consequence of the following observation: every fixed quad Q has 4 points and there are precisely 4 pairs hq;Qi with p P Q.

We can count P in another way: by fixing first a point and then counting the pairs containing it. For a point p in the cell of a dangling node, there are 2 quads containing p. If p is a point in the cell of an articulation, there are 4 quads containing p. If p is a point in an arc on the boundary of a Voronoi region of a joint (not an articulation), we get that p is in 4 quads if p is not an extremity of the arc, or p is in 2d(p) quads if it is an extremity of an arc. Here d(p) denotes the number of Voronoi regions that have p in their common boundary.

The latter paragraph can be summarized as

P = ∑ vPL S f PE v 2x v f + ∑ vPM S f PE v 4x v f + ∑ vPL S f PE v 4(x v f 1) + Θ S ; ( 17
)
where Θ S is the number of pairs containing points that are extremities of the arcs in the Voronoi diagrams. Notice that Θ S is a constant for a fixed skeleton S, thus (17) can be written as

P = ∑ vPL S f PE v 2x v f + ∑ vPM S f PE v 4x v f + ∑ vPL S f PE v 4x v f + Φ S ; ( 18
)
where Φ S = Θ S ∑ vPL S 4jE v j is a constant for a fixed skeleton S.

From equations ( 16) and ( 18) we conclude that

∑ vP(V S L S ) ∑ f PE v 2x v f + ∑ vPL S ∑ f PE v x v f ( 19
)
is an objective function that minimizes the total number of quads in a scaffold.

Integer Linear Programming models

We can now state the Integer Linear Programs (IPs) that compute subdivision for the arcs such that the scaffold have a minimal number of quads. For all the IPs the optimality criterion is to minimize the objective function [START_REF] Makhorin | GNU Linear Programming Kit -Reference Manual[END_REF]. In 3), ( 4), ( 9) (19) Regular Symmetric ( 15), ( 4), ( 9) [START_REF] Makhorin | GNU Linear Programming Kit -Reference Manual[END_REF] Table 1: Integer Linear Program models.

Theorems 3.4, 3.6 and 3.7 prove the feasibility of the standard, symmetric, and regular symmetric models respectively. As discussed in Section 3.4, the feasibility of the regular IP model is a consequence of the feasibility of the regular symmetric model.

A final observation is that once we can guarantee feasibility, and since all the variables in ( 19) are positive integers, then an optimal minimal solution always exists. The optimal solution can be computed with a Mixed-Integer Linear Program Solver. In particular with the branch-and-cut [START_REF] Chen | Applied Integer Programming[END_REF]Chapter 12] implementation of GLPK [START_REF] Makhorin | GNU Linear Programming Kit -Reference Manual[END_REF].

Algorithms

In this section we provide practical algorithms for the implementation of our method. Although our main contribution is in the theoretical foundation and the proofs of the existence of scaffolds, the algorithms in this section show that our method can be readily implemented to compute scaffolds for any skeleton.

The general steps for an implementation are described in Algorithm 1. Sub-algorithms 2, 3 and 4 deal with the details of the construction of compatible cells from the output of the standard IP and the definition of the bijections between cells (linking process).

Input: The set of nodes V S and edges E S representing the skeleton.

Output:

The quads that represent a scaffold. Algorithm 1 follows the general three steps we introduced in Section 1. In Step 1 we use spherical Voronoi diagrams for the partition of the sphere at the joints (lines 1-4). We mostly work with their duals, Delaunay triangulations, which are equivalent to the the convex hulls [START_REF] Brown | Geometric Transforms for Fast Geometric Algorithms[END_REF][START_REF] Grima | Computational Geometry on Surfaces[END_REF]. The discretization of the regions, Step 2, is divided into sub-algorithms.

1 foreach node v P V S do 2 A v 2 feS v : e P E S ; e ( vg; 3 H v 2 convexHull(A v ); 4 E v 2 edgesO f (H v );
The convex hull computed in Algorithm 1 (line 3) is not always 3-dimensional. It may be a 2-dimensional convex hull (a poly- gon), a 1-dimensional convex hull (one edge), or a 0-dimensional one (a point). Those cases correspond respectively to: more than three points in general position, more than two coplanar points, two points, and only one point (Figure 10). 1-dimensional convex hulls occur around articulations, while 0-dimensional ones are associated to dangling nodes. In our implementation the convex hull is computed by means of the QHull library [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF]. Algorithm 2 deals with the standard IP whose optimal solutions give compatible cells. It sets up and solves the IP using the convex hull representation of the regions. Simple modifications can be done for the other three variants of the scaffold. We follow ( 5). Thus we chose to have at least 4 points on each cell to guarantee at least quadrangular cross sections as in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF][START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF][START_REF] Panotopoulou | Scaffolding a skeleton[END_REF][START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF][START_REF] Yao | Compatible quadrangulation by sketching[END_REF]. We set the threshold for long arcs at 5π 6 (i.e δ = π 6 ), above which at least two subdivisions must be done. In our implementation we solve the IP with the branch-and-cut [START_REF] Chen | Applied Integer Programming[END_REF]Chapter 12] implementation in GLPK [START_REF] Makhorin | GNU Linear Programming Kit -Reference Manual[END_REF].

Input: The set of nodes V S and edges E S representing the skeleton, along with E v for each v P V S .

Output: The values x v f representing the number of subdivisions for each arc that gives compatible cells.

x a g = ∑ hPE b h((eS v ) x b h ; 12 Define objective function ∑ vP(V S L S ) ∑ f PE v 2x v f + ∑ vPL S ∑ f PE v x v f for IP;
13 Solve IP by minimizing the objective function. Algorithm 3 computes the cells from the subdivision numbers and the convex hull representation of the Voronoi diagram, and such it must handle the different cases of the convex hulls (see Figure 10).

The heuristic for the linking process (Step 3) is described in Algorithm 4, from which the resulting quads are defined. It defines the Input: Nodes V S and edges E S of the skeleton, E v ; H v and the subdivision numbers x v f .

Output: A list C of compatible cells for the scaffold. bijections (links) by minimizing the total length of the links (same heuristic used in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]).

Input: The set of compatible cells C = fC v e j v P V S ; e P E S g.

Output:

The bijections φ e between linked cells.

1 foreach edge e = ab P E S do 2

Order points in the cells C a e and C b e according to the angle around the edge e; 3 Let C a e = hp 0 ; p 2 ; : : : ; p n i and C a e = hq 0 ; q 2 ; : : : ; q n i; 4 D min 2 ∑ n i=0 kp i q i k; k 2 0; To better reflect the geometry of the model, and depending on the application at hand, the position of points in the discretization of each arc can be modified by a Laplacian smoothing [START_REF] Botsch | Polygon Mesh Processing[END_REF] or another heuristic. A global optimization could be further applied to take into account the twisting of the quads around a chain of articulation nodes (Figure 11). The challenge is in devising such a "better" positioning of points that also respects the symmetries of the skeleton. An intrinsic solution using our method is to subdivide each arc twice (or more) as shown in Figure 11.

An example of scaffold with different radii for the spheres at the joints is shown in Figure 12. Note that the skeleton is symmetric as well as the scaffold, hence the values of the radii are also symmetric. Symmetric and regular scaffolds Regular scaffolds can be obtained with the algorithms described above, a simple variation of Algorithm 2 is needed in order to account for the restrictions given in [START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF]. Similarly, for symmetric scaffolds we need to add the extra restrictions in ( 9), taking special care with the definition of links (Algorithm 4). Dangling nodes or articulations fixed by the skeleton symmetries might present issues if the linking is naively done.

In our implementation we modify Algorithm 4 such that the links are propagated from one edge to its symmetric ones. The cells of dangling nodes are computed by projecting the linked cells onto the plane perpendicular to the incident edge (of the dangling nodes). This guarantees that the links respect the symmetries.

Linear system simplification The IP in Algorithm 2 can be simplified by removing variables and equations relative to dangling nodes. The equations in a chain of articulations can be merged into one equation relating the extremal cells of the chain, removing the variables relative to intermediate articulations. In Figure 13 we illustrate the simplification cases. The general complexity of our method is bounded by the IP solver algorithm, which theoretically is NP-complete [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF]. In practice, as usual with Linear Programming, the solver behaves well. The convex hull implementation has average complexity n log n, with n the number of skeleton edges. In Table 2 we show a summary of the running time of our method for some examples. It illustrates the general relation between the total time and the time spent in solving the IP. The reported times are averages of 200 runs.

Our Python implementation is not optimal but versatile. It generates the IP problem file and call GLPK [START_REF] Makhorin | GNU Linear Programming Kit -Reference Manual[END_REF] (in particular the utility glpsol) to solve the IP. QHull [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF] is also called as an external utility and has a negligible running time. Better timings should be expected from a production-ready C++ implementation with direct linking to GLPK and QHull libraries.

6 Further simplifications of the scaffold

As we discussed in Section 3, the minimal number of points for each cell can be set to three. This decreases the number of quads in the final scaffold and might generate triangular cross-sections along some skeleton edges. In Figure 14 we recompute the example in Figure 1a with the new lower bounds, this decreases the number of quads in the standard scaffold from 192 to 144. Table 3 summarizes the number of quads and vertices per valency of some of the scaffolds in this paper, all of which are optimally computed with our implementation. It shows that the symmetric scaffolds on Figure 4 and 6 are also optimal standard scaffolds. Table 3: Summary of some scaffolds in the paper. Columns: v knumber of vertices with valency k, t quads -total number of quads.

The entries £ mean zero. Another simplification arises by noticing that any partition of the sphere with an inscribable dual gives feasible solutions for the IPs. Thus the computation of Voronoi regions can be modified by changing the minimal angle for which two triangles sharing an edge in the convex hull are considered to be coplanar. In our experience, this helps to reduce the complexity of the cells around joints and in general reduces (even more) the number of quads needed to construct a scaffold. As an example, in Figure 15 four close-to-coplanar points can be considered coplanar, yielding compatible cells with less points (hence less quads). On the other hand, due to roundoff errors coplanar points may appear as not coplanar. This issue provokes that very small arcs appear on the boundary of Voronoi regions, consequently some quads may look like triangles at a large scale (Figure 15b).

Application to polygonization

In this paper we focus on the feasibility of the scaffold construction, including the symmetric and regular variants. The scaffold is in most applications an intermediate step [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF][START_REF] Baerentzen | Interactive shape modeling using a skeleton-mesh co-representation[END_REF][START_REF] Blidia | G1-smooth splines on quad meshes with 4-split macro-patch elements[END_REF][START_REF] Ji | Bmesh: A modeling system for base meshes of 3d articulated shapes[END_REF][START_REF] Karčiauskas | Curvature continuous bi-4 constructions for scaffold-and sphere-like surfaces[END_REF][START_REF] Srinivasan | Solidifying wireframes[END_REF][START_REF] Usai | Extraction of the quad layout of a triangle mesh guided by its curve skeleton[END_REF][START_REF] Yao | Compatible quadrangulation by sketching[END_REF]. Although we present theoretical results and practical algorithms we also want to showcase an example application.

Polygonization of implicit surfaces is usually done with Marching Cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] or one of its variants [START_REF] Wenger | Isosurfaces. Geometry, Topology & Algortihms[END_REF]Chapter 2]. For computer graphics applications quad-dominant meshes are more desirable and this requires a re-meshing of the standard triangular meshes obtained with the Marching Cubes algorithms [9, Section 6.6]. For articulated models, it is a good improvement if the quad mesh follows the structure of the skeleton. A scaffolding technique can be used to bypass the re-meshing process and get directly a quad-dominant mesh that follows the skeleton.

Figure 16 shows a mesh computed with our scaffolding technique. The implicit surface (Figure 16a) is a convolution surface [START_REF] Zanni | Skeleton-based Implicit Modeling & Applications[END_REF] for which a scaffold is computed with a higher number of subdivisions for each arc: multiplying the standard solution by an integer κ > 0 that also defines the number of quads along a segment. To obtain the surface mesh we project the scaffold quads onto the surface. This is done by ray-shooting: compute the intersection with the surface of rays emanating from the skeleton in the direction given by quad vertices. A similar projection was done in [START_REF] Angelidis | Adaptive implicit modeling using subdivision curves and surfaces as skeletons[END_REF], the main difference with our approach is that the mesh we project (scaffold) is computed for the whole skeleton while in [START_REF] Angelidis | Adaptive implicit modeling using subdivision curves and surfaces as skeletons[END_REF] a fixed set of meshes is defined for each line segment. At the tips of dangling nodes, the quads are computed by creating a polar-annular region with a singular point (as done in [START_REF] Baerentzen | Converting skeletal structures to quad dominant meshes[END_REF]). The polar-annular meshes transform all valency 3 vertices on the boundary of the scaffold mesh into valency 4 vertices at the cost of having an extra high valency vertex per dangling node.

Conclusions

We presented a method to construct a coarse quad-dominant mesh around a skeleton made of line segments that can serve as intermediate step in several applications. The mesh follows the structure of the skeleton and is computed in an optimal way, minimizing the total number of quads. Our method works for any skeleton even in the presence of cycles. We modeled the problem as an Integer Linear Program and proved its feasibility. We presented variants of our method: we can generate a mesh with the same number of quads around each line segment, or a mesh with the same symmetries as the underlying skeleton. It is also possible to satisfy both conditions at the same time. The scaffold thus obtained can be used to polygonize implicit surfaces around the skeleton or as an intermediate step in other applications. Our method overcomes some limitations of, and add extra features to, previous work. The algorithmic description of each step gives a basis on which multiple implementations can be developed.

Figure 2 :

 2 Figure 2: Construction of cells: for every joint (a) take the intersection of the incident segments with the unit sphere (b), compute the Voronoi diagram (c), subdivide the arcs in the boundaries of Voronoi regions (d) and take the ordered set of points (polyline) as representation of the cells (e).

Figure 3 :

 3 Figure 3: Degenerate cases arise when allowing for only one subdivision per arc, with at least 3 (a), or 4 (b), points on each cell. In both cases long pole-to-pole arcs are discretized as single segments going through the joint.

  Ve = ab P E S :

Figure 4 :

 4 Figure 4: A symmetric scaffold is expected for a symmetric skeleton. In the picture the skeleton is symmetric through the plane π. (a) The cells of one standard scaffold does not respect the symmetry, links are not symmetric either. (b) The cells of a symmetric scaffold respect the symmetry.

T

  is a skeleton symmetry, thus T (e) P E S is the edge connecting T (a) and T (b). Since xv f is a solution of the system given by (3), taking the equation for the edge T (e) in (3), it follows that ∑ kPE T (a) k((S T (a) T(e))

  Proof. By Theorem 3.4 we have that there is a solution xv f of (3) satisfying (4). Lemma 3.5 then gives a solution xv f to (3) satisfying (9). Since xv f ! xv f , we have that (4) is trivially satisfied by xv f .

Figure 6 :

 6 Figure 6: A three-fold rotation symmetry (C 3 ). (a) Skeleton. (b) Standard scaffold. (c) Symmetric scaffold.

Figure 7 :

 7 Figure 7: Complex closed skeletons (a), standard scaffolds (b), and regular scaffolds (b)

Theorem 3 . 7 .Figure 8 :

 378 Figure 8: Symmetry and regularity are independent properties. (a) Regular asymmetric scaffold: the highlighted cell is asymmetric with respect to the reflection symmetry. (b) Symmetric irregular scaffold: one edge has a quadrangular cross profile while the rest are pentagonal.

Figure 9 :

 9 Figure 9: With symmetric requirements a scaffold can be greatly improved. (a) The computation of the scaffold in Figure 1a without symmetric requirements. (b)The "correct" symmetric scaffold respecting all the symmetries of the skeleton (same as Figure1a).

5 10 for i = 0 to n do 11 OutputAlgorithm 1 :

 10111 Define and solve the Linear Program for the subdivision of arcs; / * Alg. 2 * / 6 Construct compatible cells; / * Alg. 3 * / 7 Define bijections of linked cells; / * Alg. 4 * / 8 foreach edge e = ab P E S do 9 Let C a e = hp 0 ; p 2 ; : : : ; p n i; quad hp i ; φ e (p i ); φ e (p i+1 ); p i+1 i / * i + 1 is taken mod n * / General algorithm for constructing a scaffold.

Figure 10 :

 10 Figure 10: Cases for the convex hull of A v . Edges of E v are shown in black. (a) 3-dimensional: more than three points in general position. (b) 2-dimensional: more than two coplanar points. (c) 1dimensional: two points. (d) 0-dimensional: one point.

Algorithm 2 :

 2 Compute subdivisions for compatible cells.

H v ; 14 foreach 24 Add the points to C v e ; 25 Add C v e to C ; 26 ifAlgorithm 3 :

 142425263 node n e in H v do / * n e = e S v for e ( v * / 15 C v e 2 emptyList(); 16 Let f ; g be the two edges incidents to n e ; 17 Compute x v f + 1 points in the arc from v + N 1 to v + N 2 going perpendicularly to f on S v ; 18 Compute x v g + 1 points in the arc from v + N 1 to v + N 2 going perpendicularly to g on S v ; 19 Add the points to C v e avoiding repetitions; 20 Add C v e to C ; 21 if H v is 1-dimensional then / * only one edge and 2 nodes * / 22 Let f be the unique edge in E v ; 23 Compute x v f points in the circle with center v and perpendicular to f on S v ; H v is 0-dimensional then / * only one node * / 27 Compute x v ν v points in the great circle with center v and perpendicular to e on S v ; Construct compatible cells.

5 for j = 1 to n 1 do 6 D 2 ∑ n i=0 kp i q i+ j k; 7 if D < D min then 8 D min 2 D; k 2 j; 9 Algorithm 4 :

 56278294 Define bijection the φ e as p i U 3 q i+k ; / * i + k is taken mod n * / Construct bijections between linked cells.

Figure 11 :

 11 Figure 11: In a chain of articulations the scaffold might have a twisting behavior: (a) and (b) -different views of the same model; top row -the skeleton; middle row -the standard scaffold; and bottom row -the scaffold with extra subdivisions. Note how in the bottom row the twisting behavior of the quads is diminished.

Figure 12 :

 12 Figure 12: A scaffold with different radii for the spheres at the joints: (a) the original scaffold, (b) the scaffold with varying radii.

Figure 13 :

 13 Figure 13: Simplification of the equations in the IP.

Table 2 :

 2 Fig. 9a 133 14 119 Fig. 9b 142 19 123 Fig. 6b 55 12 43 Fig. 6c 55 12 43 Fig. 7b-top 16 4 12 Fig. 7c-top 24 4 20 Fig. 7b-mid 29 6 23 Fig. 7c-mid 37 7 30 Fig. 7b-bot 44 7 37 Fig. 7c-bot 51 6 45

  Fig. 14 * * 96 144 Fig. 6b 18 33 14 63 Fig. 1a * 48 96 192 Fig. 6c 18 33 14 63 Fig. 1b * 120 * 120 Fig. 16b 48 66 12 108 Fig. 1c * 72 * 72 Fig. 16c 96 150 12 216 Fig. 1d * 72 * 72 Fig. 16d 192 318 12 432 Fig. 4a-b 24 8 4 24 Fig. 16e 384 654 12 864

Figure 14 :

 14 Figure 14: Example of a scaffold constructed with 3 as minimal number of points on each cell. Compare with Figure 1a.

Figure 15 :

 15 Figure 15: Close-to-coplanar points can be treated as coplanar as a way to improve the scaffold. In the picture the four intersection points in the joint are not coplanar. In (a) they are treated as coplanar points, in (b) the exact Voronoi regions (and cells) are computed.

Table 1

 1

			Minimize
	Scaffold	Constraints function
	Standard	(3), (4)	(19)
	Regular	(15), (4)	(19)
	Symmetric	(

we show each model.

1

  Initialize the linear program IP;2 foreach node v P V S , and edge f P E v do

	3	Add integer variable x v f to IP;
	4 5	if the arc associated to x v f has length < 5π 6 ε v then Add restriction x v f ! 1 to IP;
	6 7	else	Add restriction x v f ! 2 to IP;
	10 foreach edge e = ab P E S do	f PE v f ((eS v )	x v f ! 4;
	11	Add following restriction to IP	∑
				gPE a
				g((eS v )

8 foreach pair (v; e) with v P V S ; e P E S and e ( v do 9 Add following restriction to IP ∑

  e in H v do / * n e = e S v for e ( v * / Compute the outward pointing unit normals N 1 ; N 2 of F 1 ; F 2 respectively; Compute x v f + 1 points in the arc from v + N 1 to v + N 2 going perpendicularly to f on S v ;

	1 C 2 emptyList(); 2 foreach v P V S do 3 if H v is 3-dimensional then
	/ * at least 4 not coplanar nodes	* /
	4 foreach node n 5 C v e 2 emptyList(); 6 foreach f ( (e S v ) do 7 Let F 1 ; F 2 be the two faces in H v that have common
	boundary f ;

8 9 10 Add the points to C v e avoiding repetitions; 11 Add C v e to C ; 12 if H v is 2-dimensional then / * at least 3 nodes, all coplanar * / 13 Let N 1 ; N 2 be the two unit normals of the plane supporting

Acknowledgments: This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675789.