Steven J Cox

Antoine Henrot

Eliciting Harmonics on Strings

Keywords: Point-wise damping, spectral abscissa, Riesz basis

One may produce the qth harmonic of a string of length π by applying the 'correct touch' at the node π/q during a simultaneous pluck or bow. This notion was made precise by a model of Bamberger, Rauch and Taylor. Their 'touch' is a damper of magnitude b concentrated at π/q. The 'correct touch' is that b for which the modes, that do not vanish at π/q, are maximally damped. We here examine the associated spectral problem. We find the spectrum to be periodic and determined by a polynomial of degree q -1. We establish lower and upper bounds on the spectral abscissa and show that the set of associated root vectors constitutes a Riesz basis and so identify 'correct touch' with the b that minimizes the spectral abscissa.

Introduction

We are interested in one's capacity to elicit harmonic overtones, or simply 'harmonics,' from stringed instruments. Such questions have long vexed theorists, though not practitioners, and were key points of friction in the development of mechanics and analysis in the period between Newton and Cauchy. Although Truesdell [START_REF] Truesdell | The Rational Mechanics of Flexible or Elastic Bodies, 1638-1788[END_REF], Cannon and Dostrovsky [START_REF] Cannon | The Evolution of Dynamics, Vibration Theory from 1687 to 1742[END_REF] and Christensen [START_REF] Christensen | Rameau and Musical Thought in the Enlightenment[END_REF] cover this period in fascinating detail, their focus on the intellectual developments of D'Alembert, D. Bernoulli, Euler, Lagrange and Rameau, leaves a number of practical innovations, notably harmonics, in the dark. As the associated mathematical question was in fact well set by 1692 and yet its mathematical analysis did not commence until 1982 we open with an outline of the full history.

Though overtones, or upper partials, had been heard through the ages they were not seen until late in the 17th century. More precisely, though Descartes, Mersenne and Galileo each remark on the sound of overtones and often ascribe them to sympathetic vibrations with neighboring bodies, Wallis [START_REF] Wallis | Concerning a new musical discovery[END_REF] in 1677 is the first to report on the means by which one may see the nodes of a vibrating string. He uses one string to excite a higher mode in a second, on which are suspended paper riders, and notes that, when the first is tuned an octave above the second, and the first is struck, "the two halves of this other [second] will both tremble; but not the middle point." He proceeds to find the nodes of what we call the third and fourth modes by tuning the first string a twelfth and double octave, respectively, above the second and observes that these nodes divide the string in equal, or aliquot, parts. Roberts [START_REF] Roberts | A discourse concerning the musical notes of the Trumpet, and trumpetmarine, and of the defects of the same[END_REF] in 1692 was the first to observe that, by softly pressing at a node while bowing elsewhere, higher modes could be elicited on a single string. He developed these ideas in his explication of the trumpet marine, or nun's fiddle, a curious variant of the ancient monochord on which the only pleasant tones are harmonics. Roberts writes "Now in the Trumpet-Marine you do not stop close as in other instruments, but touch the String gently with your thumb, whereby there is mutual concurrence of the upper and lower part of the String to produce the sound" and observes "that the Trumpet Marine will yield no musical sound but when the stop makes the upper part of the string an aliquot of the remainder, and consequently of the whole." Roberts concludes with a lovely figure that includes both a clear indication of the first three theoretical mode shapes and a practical diagram of the 15 nodes at which the player ought to gently thumb.

After Roberts the theory and practice of harmonics regrettably diverge over the course of the 18th century. On the practical side we have Mondonville's opus 4, Les Sons Harmoniques [START_REF] Banat | Masters of the Violin, v. 5, Sonatas for the Violin, Jean-Joseph Cassanéa de Mondonville[END_REF], of 1735, the first pieces scored for violin harmonics, complete with a beautiful diagram identifying notes with nodes. On the theoretical side, the mode shapes of Roberts were not to reappear until the 1747 study by Daniel Bernoulli [START_REF] Bernoulli | Réflexions et éclaircissemens sur les nouvelles vibration des cordes exposées dans les mémoires de 1747 & 1748[END_REF], a work that stemmed not from Roberts' practical considerations but from the theoretical work of Taylor [27]. In struggling to understand how a single string could simultaneously support so many modes, D. Bernoulli, Euler and Lagrange, see [START_REF] Rayleigh | [END_REF], [START_REF] Marden | Geometry of Polynomials[END_REF], and [START_REF] Cannon | The Evolution of Dynamics, Vibration Theory from 1687 to 1742[END_REF], brought forth trigonometric series but rejected these in favor of Euler and D'Alembert's general wave solution and vague reference to sympathetic resonance. None of these four men were aware of the simple mechanical means, of Roberts and Mondonville, for eliciting these individual modes. This is remarkable given the presence of J-P. Rameau, a composer and music theorist, hailed as the 'Newton of Harmony,' (see [6, p. 8]), who not only entered into substantive individual correspondence with D. Bernoulli, Euler, Lagrange and D'Alembert but who was also surely aware of both the strange flageolet tones of the trumpet marine and the works of Mondonville. Regarding Rameau's influence on Bernoulli et al. see Truesdell [START_REF] Truesdell | The Rational Mechanics of Flexible or Elastic Bodies, 1638-1788[END_REF] and especially Christensen [START_REF] Christensen | Rameau and Musical Thought in the Enlightenment[END_REF]. Galpin [START_REF] Galpin | Monsieur Prin and his trumpet marine[END_REF] informs us that Rameau's tenure with the orchestra of the Académie du Concert of Lyon overlapped with that of J.-B. Prin, the established master of the trumpet marine. Rameau, in the preface to his Pièces de Clavecin en concert, 1741, speaks of the recent success of such ensemble pieces (most likely Mondonville's opus 5, see Girdlestone [13, p. 43]) and names one of his own pieces after A.-J. Boucon, his student and the wife of Mondonville. And yet, neither Prin's nor Mondonville's harmonics figure anywhere in the vast theoretical writings of Rameau. We have not been able to discern whether this lack was intentional or whether Rameau himself failed to recognize the theoretical importance of their work. Its effect however was to postpone for 100 years any mention of harmonics and so offer support to Miller's dark verdict, [24, p. 37], that "A careful search fails to reveal any major contribution to the science of sound which arose in the 18th century."

It is not until 1807 in the work of Thomas Young [33, lect. 32] that we find a lucid synthesis of Roberts and Mondonville: "if a long chord be initially divided into any number of such equal portions, its parts will continue to vibrate in the same manner as if they were separate chords; the points of division only remaining always at rest. Such subordinate sounds are called harmonics: they are often produced in violins by lightly touching one of the points of division with the finger, when the bow is applied, and in all such cases it may be shown, by putting small feathers or pieces of paper on the string, that the remaining points of division are also quiescent, while the intervening points are in motion." Although Young's observations were echoed by Helmholtz, 1862, [15, chap. 4] and illustrated by Tyndall, 1875, [29, §3.6], it appears that Rayleigh, 1877, [25, §134] in noting that "when but a single point of the string is retarded by friction there are no normal co-ordinates properly so called," was the first to have attempted a mathematical analysis of harmonics. After Rayleigh the problem remained untouched until the beautiful, but little known, 1982 work of Bamberger, Rauch and Taylor [START_REF] Bamberger | A model for harmonics on stringed instruments[END_REF]. (Coincidentally, the problem of pointwise damping of suspension bridge cables was also first broached by Kovács [START_REF] Kovács | Zur Frage der Seilschwingungen und der Seildämpfung[END_REF] in 1982). Bamberger et al. followed Rayleigh in retarding by friction, b, at but a single point, a, the motion of a fixed string of length π with transverse displacement u. More precisely, they considered

u tt -u xx + bδ(x -a)u t = 0, u(0, t) = u(π, t) = 0,
with initial data in the energy space X = H 1 0 (0, π) × L 2 (0, π) endowed with the inner product

[f, g], [u, v] = π 0 {f u + gv} dx. (1.1)
They showed that if the damping is applied at a rational multiple of π then it only damps those modes not vanishing there. They then equated "correct touch" with the value of b that best damps these remaining modes. In the case a = π/2 they showed that b = 2 was the correct touch. We here attack the general case by investigating the spectra of the associated wave operator

A(a, b) = 0 I ∂ xx 0 x = a D(A(a, b)) = {[f, g] ∈ (H 1 0 (0, π) ∩ H 2 (0, a) ∩ H 2 (a, π)) × H 1 0 (0, π) : f (a + ) -f (a -) = bg(a)}
Bamberger et al. argue that A(a, b) has a compact resolvent and so discrete spectrum and that when a = pπ/q, where p and q are relatively prime integers, it has the imaginary eigenvalues

λ 1,n = inq, n ∈ Z \ 0 and associated eigenvectors V 1,n = sin(nqx)[1/(nq) i]. (1.2)
We write σ(pπ/q, b) for the spectrum of A(pπ/q, b) and

σ H (pπ/q, b) = {λ 1,n } (1.3)
for the harmonic part of the spectrum. It is shown in [START_REF] Bamberger | A model for harmonics on stringed instruments[END_REF] that the nonharmonic spectrum, σ(pπ/q, b) \ σ H (pπ/q, b) lies in the strict left half plane and that the semigroup e tA(pπ/q,b) enjoys exponential decay on the orthogonal complement of the harmonic subspace H ≡ span{V 1,n } More precisely, denoting the orthogonal complement of H by H ⊥ , Bamberger et al. establish the existence of positive constants, C 1 and C 2 , depending only on p/q and b, such that e tA(pπ/q,b)

U X ≤ C 1 e -C 2 t U X ∀ U ∈ H ⊥ , t ≥ 0. (1.4)
In this work we shall show the greatest C 2 for which (1.4) holds is the magnitude of the nonharmonic spectral abscissa of A(pπ/q, b),

µ(pπ/q, b) ≡ max{ λ : λ ∈ σ(pπ/q, b) \ σ H (pπ/q, b)} (1.5)
and, along the way, establish bounds on b → µ(a, b) and produce explicit minimizers when a = π/3 and a = π/4. More precisely, in §2 we delineate the structure of the spectrum, taking care to identify multiple eigenvalues and to specify the behavior of the spectrum for b near 0, 2 and ∞. In §3 we establish upper and lower bounds on the nonharmonic spectral abscissa and arrive at the correct touch for a = π/3 and a = π/4. In §4 we prove that the associated root vectors comprise a Riesz basis and so identify µ with the best rate of decay. We close in §5 with a return to dimensions and a discussion of the true magnitude of b and the role played by an additional elastic term at a. There exist many works, e.g., Cox and Zuazua [START_REF] Cox | The rate at which energy decays in a damped string[END_REF], [START_REF] Cox | The rate at which energy decays in a string damped at one end[END_REF], devoted to the optimal damping of strings, and a sizable mathematical literature, e.g., Liu [START_REF] Liu | Energy decay problems in the design of a pointwise stabilizer for string vibrating systems[END_REF], Jaffard et al. [START_REF] Jaffard | Singular internal stabilization of the wave equation[END_REF], and Ammari et al. [START_REF] Ammari | Asymptotic behavior of the solutions and optimal location of the actuator for the pointwise stabilization of a string[END_REF], and engineering literature, e.g., Kovács [18] and Krenk [START_REF] Krenk | Vibrations of a taut cable with an external damper[END_REF], that in fact focus on pointwise damping. As the focus of these previous efforts was on the optimal damping of either all modes or a single mode, the harmonics problem of Bamberger et al. was naturally neglected. (2.1)

On the Disposition of the Spectrum

The majority of our results are consequences of the simple observation that

S(λ; pπ/q, b) = -(1/4) exp(λπ)P (exp(-2λπ/q); b) (2.2)
where P is the polynomial

P (w; b) = (2 -b)w q + bw p + bw q-p -(2 + b). (2.3)
In terms of its roots, {w k = |w k | exp(iθ k )} q k=1 (repeated according to their multiplicity), the eigenvalues of A(a, b) are simply

λ k,n = -q 2π {ln |w k | + i(θ k + 2πn)}, k = 1, . . . , q, n ∈ Z (2.4)
For example, if p/q = 1/2 we find P (w; b) = (2 -b)w 2 + 2bw -2 -b. Its roots, w 1 = 1 and w 2 = (b + 2)/(b -2) translate readily into the associated eigenvalues

λ 1,n = i2n, n = 0 λ 2,n = - 1 π ln b + 2 b -2 -i(H(2 -b) + 2n)
where H is the Heaviside function. It follows that

µ(π/2, b) = - 1 π ln b + 2 b -2
approaches -∞ as b → 2, in agreement with Bamberger et al.

In preparation for the general case we now examine numerical approximation of the spectra for 0 < b < ∞ and q up to 8. Our calculations are facilitated by the fact that the spectrum is completely determined by the polynomial (2.3) and the observation that S(λ; a, b) = S(λ; a, b) and S(λ + qi, pπ/q, b) = (-1) q S(λ, pπ/q, b), from which we conclude that the spectrum has period qi and that within a period the spectrum is symmetric about iq/2.

In the next 10 panels we have plotted the first q eigenvalues as a function of b as b ascends from 0 to 5 in steps of 0.01. For b < 2 we have used left-pointing arrowheads while for b > 2 we have used right-pointing arrowheads. We have, for greater legibility, limited the plot range real parts in excess of -2. We have placed a dashed vertical line at each minimal spectral abscissa. Each plot is labeled with its particular value of a and the associated (approximate) minimal spectral abscissa, µ * , and correct touch, b * . These approximate values are obtained by solving the discrete minimization problem min{µ(a, k/100); k = 0 : 500}. We shall now prove that the root trajectories indeed behave in the simple manner illustrated in Fig. 1. Perhaps the first thing that catches one's eye is the crossing of root trajectories when p and q are odd (panels (D,F,H) and (I)). We recall that though Bamberger et al. proved that all eigenvalues are geometrically simple they neglected to consider higher algebraic multiplicity. The calculation is explicit when a = π/3. For

P (w; b) = (2 -b)w 3 + bw + bw 2 -(2 + b) = (w -1)((2 -b)w 2 + 2w + (2 + b))
and the quadratic has roots at (-

1 ± √ b 2 -3)/(2 -b), that coincide at b = √ 3. In terms of λ this states that λ = (3/(2π)) ln(2 - √ 3) + 3i/2 is a double eigenvalue of A(π/3, √ 3
). It follows that b = √ 3 is in fact the correct touch when a = π/3. We now consider, for general a = pπ/q, the nature of incoming, splitting and outgoing branches along the real axis and the line (λ) = q/2. Proposition 2.1. We suppose a = pπ/q and λ = α + iβ. (2) At p = 2m, q = 2n+1, S(α+iq/2; a, b) = i(-1) n {cosh(πα)+b cosh((π -a)α) sinh(aα)} and so b = -cosh(απ)/(sinh(aα) cosh((π -a)α)). As this is strictly increasing and takes values between 2 and ∞ the claim follows.

(3) At p = 2m -1, q = 2n + 1, S(α + iq/2; a, b) = i(-1) n {cosh(πα) + b cosh(aα) sinh((πa)α)} and so

F a (α) ≡ cosh(aα) sinh((π -a)α) cosh(πα) = - 1 b . (2.5)
It is a simple matter to confirm

F a (0) = 0, -1 ≤ F a (α) ≤ 0, lim α→-∞ F a (α) → -1/2
and that F a has a unique critical point, α a , satisfying α a < -1/2. As this critical point is the global minimizer it follows that (2. We now show that this is the only possible multiple root.

Proposition 2.2. If λ is a root of S(•; a, b) then (i) its algebraic multiplicity can not exceed 2 and (ii) if it is double then λ = q/2 mod q.

Proof: We note that S(λ; a, b) = sinh(πλ) + (b/2) cosh(πλ) -(b/2) cosh((π -2a)λ) and so

S (λ; a, b) = π 2 sinh(πλ) + π 2 (b/2) cosh(πλ) -(π -2a) 2 (b/2) cosh((π -2a)λ) = π 2 S(λ; a, b) -ba(2a -π) cosh((π -2a)λ)
and so if S(λ; a, b) = S (λ; a, b) = 0 then cosh((π -2a)λ) = 0. The latter requires that

λ = i(n + 1/2)π/(π -2a) but S(λ) = 0.
For the second part we note that λ = -(q/(2π)) log(z) where z is a zero of P . We now show that every double root of P is real. If P (z) = P (z) = 0 then qP (z) -zP (z) = 0. Dividing the latter by bq and rearranging brings

tz q-p + (1 -t)z p = (b + 2)/b, t = p/q (2.6) while z 1-q P (z) = 0 brings t/z q-p + (1 -t)/z p = (b -2)/b (2.7)
If z 1 ≡ z q-p and z 2 ≡ z p then (2.6) and (2.7) require

t (z 1 ) + (1 -t) (z 2 ) = 0 and t (z 1 ) |z 1 | 2 + (1 -t) (z 2 ) |z 2 | 2 = 0
and so either z 1 and z 2 are real or |z 1 | = |z 2 |. The latter is ruled out by the fact that |z| > 1. As z q-p and z p are both real and q -p and p are coprime it follows that z is real. This gives the desired result, since z = exp(-2λπ/q). We now give a proposition which can be viewed as a "separation" result for the different branches of the spectrum. Proposition 2.3. If a = pπ/q and b > 0 and b = 2 then the spectrum does not hit the horizontal lines λ = α + iβ where

β = k 2 , k = 1, . . . , 2q -1, k = q β = k q 2p , k = 1, . . . , 2p -1, k = p β = k q 2(q -p) , k = 1, . . . , 2(q -p) -1, k = q -p (2.8)
Proof: We claim that if z is a root of (2.3) such that either z q or z p or z q-p is real, then necessarily z is real. This will give the result since z real means that the corresponding eigenvalue λ satisfies λ = 0 or λ = q/2 or λ = q while the values of λ appearing in (2.8) correspond to z q or z p or z q-p real (and z not real) and are therefore forbidden. Now, if z p is real then

P (z; b) = z q (2 -b) + b z p + bz p -(b + 2) = 0
immediately implies that z q is real. As p and q are coprime it follows that z real. Similarly, if z q-p is real then

P (z; b) = (2 -b)z q-p z p + bz p + bz q-p -(2 + b) = 0
implies that z p is real and so z is real. Next, if z q = t is real and z p = |z| p e iθ then

P (z; b) = b|z| p e iθ + bt |z| p e -iθ + (2 -b)t -(b + 2) = 0
implies that |z| p e iθ + t|z| -p e -iθ is real. Now |z| > 1 implies that θ = 0, i.e., z p real, and so z itself is real.

Recalling fig. 1 we set r = max(p, q-p) and distinguish the exit trajectories (where the magnitude of the real part is an unbounded increasing function of 0 < b < 2) and return trajectories (where the magnitude of the real part is an unbounded decreasing function of 2 < b < ∞) from the finite loops (bounded trajectory). Note that q -r = min(p, q -p) and that q -2r < 0.

For the behavior of the branches as b → +∞, we easily see (for example by writing the polynomial equation (1-z r )(1 -z q-r ) = 2(1 -z q )/b), that the limit of all the branches are the points on the imaginary axis given by λ = im q r m = 0, . . . , r -1 λ = im q q -r m = 0, . . . , q -r -1

We will see below that the first values are the end points of the finite loops, while the other ones are the end points of return trajectories. Now, at b = 2 we have only r roots so q -r roots must escape as b ↑ 2. To follow these we write

P (z; b) = z r {(2 -b)z q-r + b + bz q-2r -(b + 2)z -r } = 0 and deduce when |z| → +∞ that z q-r ∼ b b-2 . When b ↑ 2, it follows that z ∼ b 2 -b 1/(q-r)
exp(iπ(1 + 2m)/(q -r)), m = 0, 1, . . . , q -r -1 .

Therefore, the possible imaginary parts of the exit trajectories (b ↑ 2) are β = q 2(q -r)

(1 + 2m) m = 0, 1, . . . , q -r -1 .

(2.9)

Similarly, the return trajectories (b ↓ 2) emanate from points with imaginary part

β = m q q -r m = 0, 1, . . . , q -r -1.
We now combine (2.9) with Proposition 2.3 to locate the starting points of the exit trajectories. Indeed, the first statement of (2.8) shows that there exists horizontal lines that the spectrum cannot cross every i/2. Therefore, the only possible starting point for the trajectory that exits at iβ, where β = q(1 + 2m)/(2(q -r)), is i round(β) where round returns the closest integer. This is well defined except when β is exactly half of an odd integer. Since q -r is coprime with q, this can only happen if q -r divides 2m + 1. But, since 2m + 1 < 2(q -r), the only possible case is 2m + 1 = q -r. In this case, Proposition 2.2(3) instructs us that the branch is exactly the half-line β = q/2 and it comes from both (q + 1)i/2 and (q -1)i/2. Without loss we identify the lower point, (q -1)i/2, with starting point of the loop that ends at iq/2 after splitting, and identify the upper point, (q + 1)i/2, as the starting point of the trajectory that exits the complex plane at β = q/2 after splitting. We have now proven Proposition 2.4. We set r ≡ max{p, q -p} and denote by round ± the extension of round for which round -(1/2) = 0 and round + (1/2) = 1.

(1) There exist r - These are smooth where they avoid the multiple eigenvalue. (2) There exist q -r exit trajectories,

1 proper loops, m (b), 0 ≤ b ≤ ∞, m (
ε j (b), 0 ≤ b ≤ 2, ε j (2 -) = q(1 + 2j) 2(q -r) and ε j (0) = i round + q(1 + 2j) 2(q -r) , j = 0, . . . , q -r -1
These are smooth where they avoid the multiple eigenvalue. (3) There exist q -r smooth return trajectories,

ρ j (b), 2 ≤ b ≤ ∞ ρ j (2 + ) = jq q -r
and ρ j (0) = i round jq q -r , j = 0, . . . , q -r -1

It follows that the nonharmonic spectral abscissa, b → µ(pπ/q, b) of (1.5), is continuous and bounded below on [0, ∞] and µ(pπ/q, ∞) = 0. Hence, there exists a correct touch, i.e., Corollary 2.5. For integer coprime p and q the function b → µ(pπ/q, b) attains its minimum over [0, ∞] at a finite value b * (p, q).

Perusal of the root trajectories in Fig. 1 suggests that b * occurs at a turning point (vertical tangent) of the smallest loop or perhaps along the fastest return trajectory. We now proceed to bound these loops.

Bounds on µ and b *

We apply the determinant methods of Schur and Cohn to bound the moduli of the roots of

R(w; b) = P (w; b)/(w -1) = (2 + b) p-1 j=0 w j + 2 q-p-1 j=p w j + (2 -b) q-1 j=q-p w j
The principal tool is the Schur-Cohn Criterion Proposition 3.1. (Marden [START_REF] Marden | Geometry of Polynomials[END_REF]Thm. 43,[START_REF] Ammari | Asymptotic behavior of the solutions and optimal location of the actuator for the pointwise stabilization of a string[END_REF]) The number of roots of w → R(w; b) in the ball of radius c is the number of sign changes in a sequence of q determinants,

1, ∆ 1 (b, c), ∆ 2 (b, c), ..., ∆ q-1 (b, c) where ∆ k (b, c) = det L k U k U T k L T k and L k (U k )
is the k-by-k lower (upper) triangular Toeplitz matrix composed from the lowest (highest) k coefficients of w → R(cw; b).

Starting at k = 1 we find L 1 = 2 + b and U 1 = (2 -b)c q-1 and so

∆ 1 (b, c) = (2 + b) 2 -(2 -b) 2 c 2q-2 .
In order to simplify the analysis of the higher determinants we shall assume, throughout the remainder of this section, that p = 1. In that case

L 2 = 2 + b 0 2c 2 + b and U 2 = (2 -b)c q-1 2c q-2 0 (2 -b)c q-1
and so

∆ 2 (b, c) = (b -2) 4 c 4q-4 -4(b -2) 2 c 2q + 2b 2 (4 -b 2 )c 2q-2 -4(b + 2) 2 c 2q-4 + (b + 2) 4 .
Now the turning points in S correspond to turning points in R which correspond to a double zero of one of the b → ∆ k (b, c). We detect such a zero by evaluating the associated discriminants. We denote by D j (c) the discriminant of b → ∆ j (b, c) and record D 1 (c) = 64c 2q-2 and D 2 (c) = 2 20 c 6q (c q-1 + 1) 4 (c q-1 -1) 4 (c q-4 + 8c 2q-4 + 18c q-2 + 8 + c q )• (c q-4 -8c 2q-4 + 18c q-2 -8 + c q )

We note D 1 never vanishes and that only the last factor in D 2 may vanish for c > 1.

Taking q = 4 we find that b → ∆ 2 (b, c) has a double root when c is a root of 7c 4 -18c 2 + 7. This has a unique root larger than 1, namely

c 4 = (1 + 2 √ 2)/ √ 7.
Transforming this back to the λ variables we predict a turn in one of the eigenvalues at λ = -(2/π) log(c 4 ) ≈ -0.2352 in agreement with Fig. 1. To find the associated b we must find the double root of b → ∆ 2 (b, c 4 ). This occurs at

b 4 = (6/5) √ 2
In order to determine whether this is the correct touch we must check the size of the remaining root of R. Writing R(w

; b 4 ) = (2 -b 4 )(w -w 1 )(w -w 2 )(w -w 2 ) we find R(0; b 4 ) = (b 4 -2)w 1 |w 2 | 2 or w 1 = 2 + b 4 (b 4 -2)c 2 4 = -3 -2 √ 2.
As this exceeds c 4 in magnitude, b 4 is indeed the correct touch. In the general case, a similar analysis allows us to get an easy lower bound for the spectral abscissa: Proposition 3.2. Let a = π/q, q ≥ 4 and α = (89+15 √ 17)/64 2.357, then the spectral abscissa µ(a, b) satisfies µ(a, b) ≥ -q 2π(q -2) log α .

(3.1)

Proof: As mentioned above, a turning point may arrive for a value of c such that E(c) ≡ -8c 2q-4 + c q + 18c q-2 + c q-4 -8 = 0 .

Following the rule of signs of Descartes [START_REF] Marden | Geometry of Polynomials[END_REF]Thm. 41,[START_REF] Banat | Masters of the Violin, v. 5, Sonatas for the Violin, Jean-Joseph Cassanéa de Mondonville[END_REF], we know that the number n of positive zeros and the number m of change of signs in the coefficient of E satisfy n ≤ m and m -n even. Since m = 2, E(0) < 0, E(1) > 0 and E(+∞) < 0, it follows that E has two positive zeros and only one of them being larger than 1. More precisely, since E(2) can be written E(2) = -8 * 2 2q-4 + 89 4 2 q-2 -8 and 8 * 2 q-2 > 89 4 for q ≥ 4, we see that the unique root of E, say c q , larger than 1 is smaller than 2. We use this fact to improve the bound on c q . Since

c 2 q + 18 + 1/c 2 q < 2 2 + 18 + 1/2 2 = 89/4 we get 0 = E(c q ) < -8c 2q-4 q + 89 4 c q-2 q -8
and the bound c q-2 q < α follows immediately. As the entire loop stays inside the disk of radius c q the bound (3.1) stems from the link between the modulus of w and (λ).

It is also possible to get a lower bound for the spectral abscissa in terms of b. For that purpose, let us observe that looking at ∆ 1 (b, c), we immediately obtain that it becomes negative when c crosses the value ((2 + b)/(2 -b)) 1/(q-1) , therefore the polynomial R has at least one root in the disk |w| < ((2 + b)/(2 -b)) 1/(q-1) . Of course, if we now consider the next determinant ∆ 2 (b, c), we can improve this upper bound.

We factor ∆ 2 (b, c) into

∆ 2 (b, c) = (b -2) 2 c 2(q-1) + 2(b -2)c q + 2(b + 2)c q-2 -(b + 2) 2 (b -2) 2 c 2(q-1) -2(b -2)c q -2(b + 2)c q-2 -(b + 2) 2 ≡ Q 1 (b, c)Q 2 (b, c)
and look for the first zero of either

Q 1 or Q 2 . We observe that Q 2 (b, c) ≤ Q 1 (b, c) while (b -2)c q + (b + 2)c q-2 ≥ 0, i.e. Q 2 (b, c) ≤ Q 1 (b, c) ⇐⇒ b ≥ 2 or b < 2 and 1 ≤ c ≤ 2 + b 2 -b 1/2 . Now, since Q 1 (b, 1) < 0 and Q 1 (b, 2 + b 2 -b 1/2 ) = (b + 2) 2 2 + b 2 -b q-3 -1 > 0
the first zero of Q 1 (and therefore of P ) is less than ((2 + b)/(2 -b)) 1/2 and, therefore, we can assume that we are in the case

Q 2 (b, c) ≤ Q 1 (b, c). At last, since Q 1 (b, 1) = -4b < 0 and Q 2 (b, 1) = -12b < 0, Q 1 (b, c) will vanish before Q 2 (b, c
), so we can restrict ourselves to the study of Q 1 (b, c).

We begin with the case b ≥ 2. Let z 0 be the positive root of

R 2 (z) ≡ (b -2) 2 z 2 + [2(b -2) + 2(b + 2)]z -(b + 2) 2 Namely, z 0 = (-2b + √ b 4 -4b 2 + 16)/(b -2) 2 for b = 2, and z 0 = 2 for b = 2. We note that z 0 > 1, set c 0 ≡ z 1/(q-2) 0 > 1 and observe that if b ≥ 2 then Q 1 (b, c 0 ) = (b -2) 2 c 2 0 z 2 0 + 2(b -2)c 2 0 z 0 + 2(b + 2)z 0 -(b + 2) 2 ≥ (b -2) 2 z 2 0 + 2(b -2)z 0 + 2(b + 2)z 0 -(b + 2) 2 = 0 .
Therefore, c 0 is an upper bound for the first zero of P when b ≥ 2.

We use the same kind of method in the case b < 2. We introduce z 1 the (unique) positive root of

R 3 (z) ≡ (b -2) 2 z 3 + 2(b -2)z 2 + 2(b + 2)z -(b + 2) 2 = 0 . Since R 3 (1) = -4b < 0 and R 3 (2/(2 -b)) = b 2 (b + 2)/(2 -b) > 0 we see that 1 < z 1 < 2/(2 -b). We set c 1 ≡ z 1/(q-2) 1 and observe that if b ≤ 2 then Q 1 (b, c 1 ) = (b -2) 2 c 2 1 z 2 1 + 2(b -2)c 2 1 z 1 + 2(b + 2)z 1 -(b + 2) 2 = c 2 1 z 1 (b -2)[(b -2)z 1 + 2] + 2(b + 2)z 1 -(b + 2) 2 ≥ z 2 1 (b -2)[(b -2)z 1 + 2] + 2(b + 2)z 1 -(b + 2) 2 = 0
where we used the fact that (b -2)z 1 + 2 > 0 and c 2 1 ≤ z 1 (for q ≥ 4). Therefore, c 1 is an upper bound for the first zero of P when b < 2. Now, using the link between the spectral abscissa µ and the moduli of the roots of R we finally obtain the following result:

Proposition 3.3. Let z 1 (b) be the (unique) positive root of R 3 (z) ≡ (b -2) 2 z 3 + 2(b -2)z 2 + 2(b + 2)z -(b + 2) 2 and z 0 (b) be the positive root of R 2 (z) ≡ (b -2) 2 z 2 + 4bz -(b + 2) 2 . Then, for q ≥ 4 µ(π/q, b) ≥ - q 2π(q -2) log |z 1 (b)| when b < 2 log |z 0 (b)| when b ≥ 2
The estimate is particularly good for large values of b as shown in the following picture drawn in the case q = 10. In order to get bounds for "the optimal touch" b * , we now need an upper bound for the spectral abscissa, at least for a particular choice of b. Since we observe numerically that b * is not so far from the value q/3, we choose this value. Proposition 3.4. If q ≥ 6 then µ(pπ/q, q/3) ≤ -3/(8πq).

Proof: We outline the key steps but suppress the many routine calculations. We recall the S of (2.2) and set φ(α, β) ≡ |S(α + iβ)| 2 . We begin by working with β fixed in the range of possible values corresponding to the different branches of the spectrum. We denote by B this set of possible values for β. We will often need to separate the cases of the closed loops and branches returning from infinity.

Step 1: The function α → φ α (α, β) is non decreasing for -1/(πq) ≤ α ≤ 0 (we just need to compute the second derivative φ αα (α, β) which is clearly non negative).

Step 2: The inequalities φ α (0, β) ≥ 0 and φ α (-1/(πq), β) + φ α (0, β) ≥ 0 hold. Actually, it is the longer and more technical part of the proof. For that, we write

1 π φ α -1 πq , β + φ α (0, β) = A + B 1 p 3 cos 2(q -p) q βπ + B 2 q -p 3 cos 2p q βπ (3.2)
with

A = q 3 + q 3 cosh 2 q -(1 + q 2 36 ) sinh 2 q - q(q -2p) 36 sinh 2(q -2p) q 2 B 1 = -1 + q 6 sinh 2p q 2 -cosh 2p q 2 B 2 = -1 + q 6 sinh 2(q -p) q 2 -cosh 2(q -p) q 2
and we get the result using classical inequalities satisfied by the functions cosh and sinh and precise study of the function of β defined in (3.2).

Step 3: We deduce from steps 1 and 2 that |φ α (α, β)| ≤ φ α (0, β) for -1/(πq) ≤ α ≤ 0 and so φ(0, β) -φ(α, β) ≤ φ α (0, β)|α| there and then, φ(α, β) (and also S) cannot vanish while

φ(0, β) -φ α (0, β)|α| > 0 ⇐⇒ |α| < φ(0, β) φ α (0, β) ≡ ψ(β)
As a consequence, the spectral abscissa µ satisfies µ ≤ -ψ(β) for each β ∈ B.

Step 4: We find a lower bound for ψ(β) = sin 2 βπ + q 2 9 sin 2 p q βπ sin 2 (q-p) q βπ 2π p 3 sin 2 q-p q βπ + q-p 3 sin 2 p q βπ . (

when β varies in the range of possible values. We now examine separately the case of a closed loop and the case of a branch coming back from infinity (since b = q/3 ≥ 2, we do not need to look at the case of a branch going to infinity). For sake of simplicity, we assume here that r = q -p (i.e. q ≥ 2p). The proof would be exactly the same in the other case.

The case of a loop finishing at imq/(q -p) corresponds to values of β belonging to the interval k, mq q-p with k ≡ round(mq/(q -p)) and m are connected by a relation like mq -k(q -p) = m 0 . Now, we use that on such an interval sin 2 ((π -a)β) ≤ sin 2 (aβ) to get ψ(β) ≥ 3 2πq

q 2 9 sin 2 ((π -a)β) + sin 2 βπ sin 2 (aβ) . (3.4)
We split the interval I = k, mq q-p into

I = k, k + 1 2 m 0 q-p ∪ k + 1 2 m 0 q-p , k + m 0 q-p = I 1 ∪I 2 . On I 1 , we have ψ(β) ≥ 3 2πq q 2 9 sin 2 q -p q βπ ≥ q 6π sin 2 q -p q k + 1 2 m 0 q -p π = q 6π sin 2 m 0 π 2q . At last, since on [0, π/12], sin 2 x ≥ 12x π sin π 12 2 , we get sin 2 m 0 π 2q ≥ sin 2 π 2q ≥ 144 π 2 π 2 4q 2 sin 2 π 12 = 18 q 2 (1 - √ 3 2 ), therefore, on I 1 ψ(β) ≥ 3(2 - √ 3) 2πq . (3.5)
On I 2 , we write:

ψ(β) ≥ 3 2πq sin 2 βπ sin 2 p q βπ ≥ 3 2πq sin 2 (kπ + 1 2 m 0 π q-p ) sin 2 pπ q (k + m 0 q-p ) , the right-hand side is equal to 3 2πq sin 2 m 0 π 2(q-p) sin 2 m 0 π (q-p) = 3 2πq 1 4 cos 2 m 0 π 2(q-p) therefore, on I 2 ψ(β) ≥ 3 8πq . (3.6)
In the case of a branch coming back from infinity, which corresponds to an interval like k, jq q-r with k ≡ round(jq/(q -r) and jq -k(q -r) = q -j 0 , we use the same idea and we get the same inequalities as (3.5) and (3.6).

As a consequence of the upper and lower bounds on the nonharmonic spectral abscissa we obtain the following bounds on the correct touch, Proposition 3.5. Let us denote by b * (p, q) the correct touch, i.e., the minimizer of b → µ(pπ/q, b). If p = 1 and q ≥ 6 then

b 1 (q) ≤ b * (p, q) ≤ b 2 (q) (3.7)
where

b 1 (q) = 2z 3 -z 2 -z + 2 + √ -7z 4 + 18z 3 -7z 2 z 3 -1 b 2 (q) = 2(z 2 -z + 1) + 2 2z(2 -z)(z -1 2 ) z 2 -1
and z = exp 3(q -2)/4q 2 .

Proof: To find these bounds, it suffices to look at the intersection of the curve giving the lower bound with the line µ = -3/(8πq). Or, equivalently, to find b such that (in the language of Proposition 3.3) z 0 (b) and z 1 (b) equal exp 3(q -2)/4q 2 . We obtain the result by plugging this value of z in the equation defining z 0 (b) and z 1 (b) respectively and solving in b.

Regarding the Root Vectors

Recall that Rayleigh was already aware of the fact that pointwise damping destroys the "normal co-ordinates" or orthonormal basis of eigenmodes enjoyed by the undamped wave operator. For our purposes however it will suffice to show that the (generalized) eigenmodes, or root vectors, are isomorphic to an orthonormal base. In other words, we show that these modes comprise a Riesz basis for the energy space, X. We begin with a specification of these modes.

In addition to the harmonic modes in (1.2) we have the nonharmonic modes

V k,n (x) = y k,n (x)[1/λ k,n 1], y k,n (x) = sinh(λ k,n (π -a)) sinh(λ k,n x) if x < a, sinh(λ k,n a) sinh(λ k,n (π -x)) if x > a, (4.1)
where k = 2, . . . , q and n runs through the integers.

If λ k,n = λ k+1,n then, in addition to its eigenvector V k,n , there exists the generalized eigenvector V k+1,n = [ỹ z] for which (A -λ k )V k+1,n = V k,n . In components this reads z = λ k,n ỹ + y k,n /λ k,n and ỹ -(λ k,n + bδ a )z = y k .

The first into the latter yields

ỹ -λ 2 k,n ỹ -λ k,n bδ a ỹ = (2 + bδ a /λ k,n )y k,n
that is, ỹ is the

H 1 0 (0, π) function satisfying ỹ -λ 2 k,n ỹ = 2 sinh(λ k,n (π -a)) sinh(λ k,n x) if x < a, sinh(λ k,n a) sinh(λ k,n (π -x)) if x > a, (4.2)
and the transmission condition

ỹ (a + ) -ỹ (a -) = λ k,n bỹ k,n (a) + by k,n (a)/λ k,n (4.3)
The solution to (4.2)-(4.3) is, piecewise,

ỹ(x) = 1 λ k,n (π -a) cosh(λ k,n (π -a)) sinh(λ k,n x) + x cosh(λ k,n x) sinh(λ k,n (π -a)) a cosh(λ k,n a) sinh(λ k,n (π -x)) + (π -x) cosh(λ k,n (π -x)) sinh(λ k,n a)
We now show that the root vectors comprise a Riesz Basis for X. In some related cases, e.g., Kergomard et al. [START_REF] Kergomard | Resonance modes in a 1-D medium with two purely resistive boundaries: calculation methdos, orthogogonality and completeness[END_REF], the underlying isomorphism can be constructed by hand. We follow here the indirect approach of Bari, e.g., Young [32,Thm. 1.9], and show that they comprise a complete Bessel sequence and are biorthogonal to a set that is also complete and Bessel. To begin, as {V 1,n } n =0 is an orthonormal sequence in X, Bessel's inequality reads

±∞ n=±1 | V 1,n , V | 2 ≤ V 2 ∀ V ∈ X
and so we need only show

Proposition 4.1. If V ∈ X then n | V k,n , V | 2 < ∞ for each k = 2, . . . , q.
Proof: We write V = [y z] and find

V k,n , V = (1/λ k,n ) π 0 y k,n y dx + π 0 y k,n z dx. (4.4)
Now, as λ k,n = λ k,0 -iqn,

1 λ k,n π 0 y k,n y dx = (-1) n(q-p) sinh(λ k,0 (π -a)) a 0 cosh((λ k,0 -inq)x)y (x) dx- (-1) np sinh(λ k,0 a) π a cosh((λ k,0 -inq)(π -x))y (x) dx (4.5)
As {nq} n is a separated sequence it follows from Young [START_REF] Young | An Introduction to Nonharmonic Fourier Series[END_REF]Thm. 4.4] that e inq is a Bessel sequence in L 2 (-π, π) and hence the sequence (4.5) is square summable in n. This same argument allows us to conclude that the second term in (4.4) is also square summable in n. Finally, given our explicit representation of the root vectors associated with double eigenvalues, the argument above may be used to show that these too constitute a Bessel sequence.

These same estimates hold for the adjoint system

A * (a, b) = 0 -I -∂ xx 0 x = a D(A * (a, b)) = {[f, g] ∈ (H 1 0 (0, π) ∩ H 2 (0, a) ∩ H 2 (a, π)) × H 1 0 (0, π) : f (a + ) -f (a -) = -bg(a)} for its root vectors are W k,n (x) = y k,n (x)[1/λ k,n -1]. It is a simple matter to confirm that this sequence is biorthogonal to {V k,n }.
Regarding completeness, as our pointwise damping is an unbounded 'perturbation' of the skew wave operator we may not argue as in Cox and Zuazua [START_REF] Cox | The rate at which energy decays in a damped string[END_REF]. Instead, following Krein and Nudelman [START_REF] Krein | On direct and inverse problems for the boundary dissipation frequencies of a nonuniform string[END_REF], Cox and Zuazua [START_REF] Cox | The rate at which energy decays in a string damped at one end[END_REF] and Ammari et al. [START_REF] Ammari | Asymptotic behavior of the solutions and optimal location of the actuator for the pointwise stabilization of a string[END_REF], we invoke the simple trace criterion of Livšic. where the ν n are repeated according to their algebraic multiplicity. Equality holds in (4.6) if and only if the root vectors of T are complete in H.

It is now a simple matter to confirm that equality holds in (4.6). We note that an alternate route to completeness has recently been established in the works of Guo and Xie [START_REF] Guo | A sufficient condition on Riesz basis with parenthesis of nonself-adjoint operator and application to a serially connected string system under joint feedbacks[END_REF] and Xu and Guo [START_REF] Xu | Riesz basis property of evolution equations in Hilbert spaces and application to a coupled string equation[END_REF].

Proposition 4.3. tr(A -1 (a, b) ) = -b(π -a)a/π. Proof: If [f g] = A -1 (a, b)[u v] then g = u and f = v + bδ a u with f (0) = f (π) = 0. Hence f (x) = x 0 (x -t)v(t) dt -(x/π) π 0 (π -t)v(t) dt -(b/π)u(a)ψ(x) where ψ(x) = (π -a)x if x < a a(π -x) if x > a Similarly, If [f g] = A -1 (a, b) * [u v] then g = -u and f = -v +bδ a u with f (0) = f (π) = 0. Hence, as above f (x) = - x 0 (x -t)v(t) dt + (x/π) π 0 (π -t)v(t) dt -(b/π)u(a)ψ(x).
Combining these calculations we find A The polynomial representation of our shooting function affords the elementary factorization S(λ; pπ/q, b) = -1 4 exp(λπ)P (exp(-2πλ/q); b)

= 1 4 exp(λπ)(b -2) q k=1 (exp(-2πλ/q) -w k ) (4.7)
The beauty of this representation is that sums of eigenvalues, as required by (4.6), appear on differentiation with respect to λ. The computation of these sums will be facilitated by

n 1 (n + α) 2 + β 2 = π 2β sinh(2πβ) cosh 2 (πβ) -cos 2 (πα) (4.8)
that follows by summing the residues of -π cot(πz)/((z + α) 2 + β). (e -2πλ/q -w k ) -π 2q e λπ(1-2/q) (b -2) q j=1 k =j (e -2πλ/q -w k )

The first term vanishes at λ = 0 and the products in the second term each vanish except when j = 1. Hence

S λ (0; a, b) = π 2q (2 -b) q k=2 (1 -w k ) = π 2q P (w) w -1 w=1 = π
as claimed. Proceeding to the second derivative

S λλ (λ; a, b) = π 2 4 e λπ (b -2) q k=1 (e -2πλ/q -w k ) - π 2 (q -1) q 2 e λπ(1-2/q) (b -2) q j=1 k =j (e -2πλ/q -w k ) + π 2 q 2 e λπ(1-4/q) (b -2) q j=1 q i=1 k ∈{i,j} (e -2πλ/q -w k ) At λ = 0 this becomes S λλ (0; a, b) = π 2 (q -1) q 2 (2 -b) q k=2 (1 -w k ) + 2π 2 q 2 (b -2) q j=2 k ∈{1,j} (1 -w k ) = 2π 2 (q -1) q - 4π 2 q q k=2 1 1 -w k = 2π 2 q q k=2 1 -w k 1 -w k - 2 1 -w k = 2π 2 q q k=2 w k + 1 w k -1 (4.9)
and so it remains only to link the w k back to their λ k,n . Well, it follows directly from (2.4) that In order to reconcile (4.10) and (4.9) we note that if w k is real then

(1/λ k,n ) = -log |w k | 2πq 1 log 2 (|w k |)/(4π 2 ) + (n + θ k /(2π)) 2 With α = θ k /(
1 -|w k | 2 |w k | 2 -2 w k + 1 = w k + 1 1 -w k while if w k is not real then 1 -|w k | 2 |w k | 2 -2 w k + 1 + 1 -|w k | 2 |w k | 2 -2 w k + 1 = w k + 1 1 -w k + w k + 1 1 -w k
As each w k is either real or occurs with its complex conjugate, we have shown, (1/λ k,n ) = tr((A -1 ) ).

S
It follows that the root vectors are complete in X and hence comprise a Riesz basis in X, and, more to the point, that the set of nonharmonic root vectors {V k,n : k = 2, . . . , q; n ∈ Z} is a Riesz basis for H ⊥ . From here it is easy to see, e.g., [START_REF] Cox | Aye there's the rub, An inquiry into how a damped string comes to rest, Six Themes on Variation[END_REF], that the nonharmonic spectral abscissa, (1.5), indeed coincides with the decay rate on H ⊥ .

Dimensions and Extensions

Though it has been convenient to work in the setting of Bamberger et al. it remains to see whether the 'correct' touch is in fact practical. To see this we start from the dimensional problem T u xx -ρu tt + γδ(x -p /q)u t = 0, 0 < x < .

(5.1)

In terms of the new space and time variables, y = πx/ and τ = (πt/ ) T /ρ we recover our initial equation u yy -u τ τ + bδ(y -pπ/q)u τ = 0, 0 < y < π where b = γ √ T ρ

We note that Cuzzucoli and Lombardo [START_REF] Cuzzucoli | A physical model of the classical guitar, including the player's touch[END_REF] consider finger damping up to 8 √ T ρ. On an absolute scale, if we take, e.g., the string studied in Cox [START_REF] Cox | Aye there's the rub, An inquiry into how a damped string comes to rest, Six Themes on Variation[END_REF], where T = 27 m kg/s 2 and ρ = 0.0015 kg/m we find that √ T ρ ≈ 0.2 kg/s corresponds to light damping.

Finally, following, e.g., Birch and Srinivasan [START_REF] Birch | Experimental determination of the viscoelastic properties of the human fingerpad[END_REF], we note that the interaction of string of finger is much more complex than our 1-parameter pointwise damper. From the 15parameter model of [START_REF] Birch | Experimental determination of the viscoelastic properties of the human fingerpad[END_REF] we may abstract an effective fingertip stiffness, k > 0, and so consider u tt -u xx + δ a (bu t + ku) = 0 (5.2)

This breaks the lovely symmetry of the shooting function, S(λ; a, b) = λ sinh(λπ) + (bλ + k) sinh(λa) sinh(λ(π -a)).

(5.3) and so also destroys the simple structure of the spectrum. In comparing figures 1 and 3 we see at once both the unfolding of the multiple eigenvalue and the genesis of the exit and return trajectories. Although the qi periodicity is broken, it appears that the full spectral abscissa is attained on one of the first q trajectories. It remains of course to determine whether the correct touch in the sense of Bamberger, Rauch and Taylor is indeed the most pleasing to the ear. If in fact the nonharmonic spectral abscissa is the right measure to minimize then Fig. 1 seems to indicate a number of improvements over current practice. More precisely, modern practice, as dictated by Zukovsky [START_REF] Zukovsky | On violin harmonics[END_REF], has the musician lightly finger at /q and yet, for odd q > 3 Fig. 1 states that fingering at 2 /q is better, i.e., µ(2 /q, b * (2, q)) < µ( /q, b * (1, q)). By this reasoning it appears that touching simultaneously at /q and 2 /q should in turn further enhance the effect.

  If A(a, b)[y z] = λ[y z] then z = λy and y must satisfy the quadratic eigenvalue problem y (x) -λ 2 y(x) -λbδ(x -a)y(x) = 0, y(0) = y(π) = 0. It is not difficult to show that y(x) = sinh(λ(π -a)) sinh(λx) if 0 < x < a sinh(λa) sinh(λ(π -x)) if a < x < π and λ is a zero of the shooting function S(λ; a, b) = sinh(λπ) + b sinh(λa) sinh(λ(π -a)).

5 ( 11 Figure 1 .

 5111 Figure 1. Eigen-trajectories as a function of touch, b, for ten choices of site, a.

( 1 )( 2 )( 3 )

 123 (0) If b < 2 there exist no real eigenvalues. If b > 2 then there exists precisely one real eigenvalue α 0 (b). It is strictly increasing with b and obeys lim b↓2 α 0 (b) = -∞ and lim b→∞ α 0 (b) = 0. If q is even and b < 2 then α 1 (b) +iq/2 is an eigenvalue and α 1 (b) is strictly decreasing with b and obeys lim b↑2 α 1 (b) = -∞ If q is odd and p is even and b > 2 then α 2 (b) + iq/2 is an eigenvalue and α 2 (b) is strictly increasing with b and obeys lim b↓2 α 2 (b) = -∞ and lim b→∞ α 2 (b) = 0. If q is odd and p is odd then there exists a unique b * < 2 and α(b * ) such that a double eigenvalue occurs at α(b * ) + iq/2. As b exceeds b * this eigenvalue splits into the pair α ± (b) + iq/2 satisfying α -(b) < α(b * ) < α + (b), lim b↑2 α -(b) = -∞ and lim b→∞ α + (b) = 0.Proof: (0) If λ = α is real then (2.1) reads b = -sinh(απ)/(sinh(aα) sinh((π -a)α)). As this is strictly increasing in α and takes values between 2 and ∞ the claim follows.(1) If q = 2n then S(α + in; p/(2n), b) = (-1) n {sinh(πα) + b cosh(aα) cosh((π -a)α)} and so b = -sinh(απ)/(cosh(aα) cosh((π -a)α)). As this is strictly decreasing in α and takes values between 0 and 2 the claim follows.

  5) has (i) no solution for b < b a ≡ -1/F a (α a ), (ii) one double root when b = b a , (iii) exactly two roots when b a < b < 2 (with one going to -∞ as b → 2, and (iv) exactly one root when b > 2.

  ∞) = imq/r, and m (0) = i round -(mq/r) m = 1, . . . , r -1

Figure 2 .

 2 Figure 2. Estimate from above (solid line) for the spectral abscissa (dashed line) in the case a = π/10.In order to get bounds for "the optimal touch" b * , we now need an upper bound for the spectral abscissa, at least for a particular choice of b. Since we observe numerically that b * is not so far from the value q/3, we choose this value.

Proposition 4 . 2 .

 42 (Krein and Langer [19, §2.5]) If H is Hilbert and T : H → H is linear and compact and T ≡ (T + T * )/2 is nonpositive and of finite trace then tr(T ) ≤ ν n ∈σ(T)

  -1 (a, b) [u v] = [-(b/π)u(a)ψ 0], and so the range of A -1 (a, b) is spanned by the unit vector e 1 ≡ [ψ 0]/ aπ(π -a). The claim follows from trA -1 (a, b) = A -1 (a, b) e 1 , e 1 .

Proposition 4 . 4 .

 44 If b ≥ 0 and b = 2 then the shooting function satisfiesS(λ; pπ/q, b) = πλ -πλ 2 q k=1 n (1/λ k,n ) + O(λ 3 )Proof: We simply differentiate (4.7) with respect to λ and make frequent use of the fact that w 1 = 1. To begin, with a = pπ/q, S λ (λ; a, b) = π 4 e λπ (b -2) q k=1

π q 1 -

 1 2π) and β = log(|w k |)/(2π) it now follows from (4.8) that n(1/λ k,n ) = -log |w k | 2πq π 2 log |w k | sinh(log |w k |) cosh 2 (log(|w k |)/2) -cos 2 (θ k /2) = -π 4q |w k | -|w k | -1 (|w k | + 2 + |w k | -1 )/4 -cos 2 (θ k /2) = |w k | 2 |w k | 2 -2 w k + 1 (4.10)

  λλ (0; a, b) = -2π q k=1 n (1/λ k,n ) as claimed. Comparing now the MacLaurin developments of S as represented in (2.1) and the above proposition, we find πλ + b(π -a)aλ 2 + O(λ 3 ) = πλ -πλ 2 λ n ∈σ(a,b) (1/λ n ) + O(λ 3 ) and so, recalling proposition 4.2, we have shown that, when b = 2, λ k,n ∈σ(a,b)

Figure 3 .

 3 Figure3. A few root paths of (5.3) when p/q = 1/5, k = 0.2, and b runs from 0 to 8.


A I Ginnis
email: ginnis@naval.ntua.gr

K V Kostas

P D Kaklis

Construction of Smooth Branching Surfaces using T-splines

Keywords: Branching surface, T-splines, Voronoi diagram, surrounding curve, discrete Fréchet distance

The request for designing or reconstructing objects from planar cross sections arises in various applications, ranging from CAD to GIS and Medical Imaging. The present work focuses on the "one-to-many" branching problem, where one of the planes can be populated with many, possibly tortuous and densely packed, contours. The proposed method combines the proximity information offered by the Euclidean Voronoi diagram with the concept of surrounding curve, introduced in [1], and T-splines technology [2] for securing a flexible and portable representation. Our algorithm delivers a single T-spline that deviates from the given contours less than a user-specified tolerance, measured via the so-called discrete Fréchet distance [3] and is C 2 everywhere except from a finite set of point-neighborhoods. Subject to minor enrichment, the algorithm is also capable to handle the "many-to-many" configuration as well as the global reconstruction problem involving contours on several planes.

Introduction

Designing or reconstructing objects from planar cross sections arises in various application areas, ranging from Computer-Aided Design (CAD) [START_REF] Pei | Constructive G 1 connection of multiple freeform pipes in arbitrary poses[END_REF][START_REF] Bentamy | Cross-sectional design with curvature constraints[END_REF] to Geographic Information Systems (GIS) [START_REF] Fujimura | Shape reconstruction from contours using isotopic deformation[END_REF][START_REF] Sirakov | A system for reconstructing and visualising three-dimensional objects[END_REF] and Medical Imaging [START_REF] Zou | Topology-constrained surface reconstruction from cross-sections[END_REF][START_REF] Schmitt | New methods for the computerassisted 3D reconstruction of neurons form confocal image stacks[END_REF][START_REF] Moschino | Joint registration and averaging of multiple 3D anatomical surface models[END_REF]. Developing generic methods for handling the user needs in the aforementioned application areas poses a variety of challenges of topological, geometrical and modeling nature, such as complying with the genus of the underlying surface, securing the required smoothness and accuracy and, finally, relying on as portable as possible surface representations.

State-of-the art algorithms have so far adopted a divide-and-conquer strategy within the cells created by the arrangement of the given cross-sections and may be coarsely categorized into a pair of families. The first one adopts surface extraction through volume reconstruction via, e.g., volume rendering, in particular marching cubes, pruning the Delaunay triangulation between adjacent cross sections [START_REF] Boissonnat | Shape reconstruction from planar cross sections[END_REF], implicit models [START_REF] Savchenko | Function representation of solids reconstructed from scattered surface points and contours[END_REF][START_REF] Turk | Shape transformation using variational implicit forms[END_REF], set-valued functions [START_REF] Cohen-Or | Guided multi-dimensional reconstruction from cross-sections[END_REF] and multi-level partition-ofunity implicit models [START_REF] Braude | Contour-based surface reconstruction using MPU implicit models[END_REF]. The second family includes interpolation approaches based on splines, in most cases low-degree (linear) piecewise interpolants (triangulations), involving global optimality criteria, such as bounding the maximum volume [START_REF] Keppel | Approximating complex surfaces by triangulation of contour lines[END_REF], minimizing the area of the surface [START_REF] Fuchs | Optimal surface reconstruction from planar contours[END_REF][START_REF] Meyers | Surfaces from contours[END_REF][START_REF] Sloan | Pessimal guesses may be optimal: A counterintuitive search result[END_REF], simplifications of the generalized Voronoi diagram [START_REF] Oliva | 3D reconstruction of complex polyhedral shapes from contours using a simplified generalized Voronoi diagram[END_REF], tiling with the aid of immersed medial axis [START_REF] Bajaj | Arbitrary topology shape reconstruction from planar cross sections[END_REF], straight-skeletons of the symmetric difference of consecutive slices [START_REF] Barequet | Nonlinear interpolation between slices[END_REF], and combining trimming of contour-surrounding curves with hole filling [START_REF] Gabrielides | G 1 -smooth branching surface construction from cross sections[END_REF]. In the same category, [START_REF] Kels | Reconstruction of 3d objects from 2d cross-sections with the 4-point subdivision scheme adapted to sets[END_REF] proposes a method for the iterative refinement of sets, corresponding to cross-sections, producing weighted averages of sets that can be used, through a subdivision scheme, for the representation of the 3D object.

The cases of non-parallel cross sections [START_REF] Boissonnat | Shape reconstruction from unorganized cross-sections[END_REF] or non-manifold surface networks that partition the space into regions with multiple labels [START_REF] Barequet | Reconstruction of multi-label domains from partial planar crosssections[END_REF] are also attracting increasing interest in pertinent literature. Finally, one should point to recent contributions which can preserve the known genus of the original object by, e.g., appropriate sampling conditions [START_REF] Amini | Geometric Tomography with Topological Guarantees[END_REF] or enriching the classical divideand-conquer strategy with global combinatorial optimization for exploring feasible topologies of tiles and computing a score that assesses the likelihood of each topology [START_REF] Zou | Topology-constrained surface reconstruction from cross-sections[END_REF].

The present work develops a method that is able to deliver smooth solutions to the "oneto-many" local reconstruction problem for cases where one of the planes (plane-1; see Figure 1) is populated with many, possibly tortuous and densely packed, contours. The paper improves drastically the method presented in [START_REF] Gabrielides | G 1 -smooth branching surface construction from cross sections[END_REF] by strengthening its topological robustness and expanding its representation portability.

The topological structure is built upon the proximity information contained in the planar Euclidean Voronoi diagram (EVd) of the contours set on plane-1. More accurately, we firstly improve the flexibility of the convex-hull based concept of the surrounding curve, introduced in [START_REF] Gabrielides | G 1 -smooth branching surface construction from cross sections[END_REF], so that, on the basis of EVd and a user-defined viewing angle ω, it is automatically decided which contours of plane-1 will be connected to parts of the single contour in the other plane (plane-0). In the sequel, we generate an EVd-based quadrilateral tiling of a 2D circular domain with holes (see Figure 2), which is homeomorphic to the sought-for surface. This tiling determines the number and the neighborhood of the surface's saddle points; see Figure 3.

Regarding the representation portability, we employ T-splines technology [START_REF] Sederberg | T-splines and TNURCCs[END_REF], which offers a flexible representation that is NURBS-compatible but free from deficiencies stemming from their non-local tensor-product nature. T-splines technology is among the recently developed NURBS generalizations that allow local refinements, such as THB-splines [START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF], LR B-splines [START_REF] Dokken | Polynomial splines over locally refined box-partitions[END_REF], Hierarchical B-splines and T-splines. In comparison with [START_REF] Gabrielides | G 1 -smooth branching surface construction from cross sections[END_REF] that delivers a G 1 -continuous multi-patch NURBS surface with trimmed patches, the present method guarantees that the output is a single T-spline that is C 2 everywhere except from a finite set of point-neighborhoods, where smoothness degrades to G 1 level. Furthermore, our numerical experience indicates that the fairness of the surfaces resulting from the proposed method is superior in comparison with those delivered from [START_REF] Gabrielides | G 1 -smooth branching surface construction from cross sections[END_REF].

Based on the "one-to-many" reconstruction process, we can handle the more general "many-tomany" case, by introducing, similarly to the approach followed in [START_REF] Park | 3-D shape reconstruction from 2-D cross-sections[END_REF], a contour on an intermediate plane and applying the "one-to-many" method twice. In our case, the introduced contour can be constructed as the intersection of the intermediate plane with the surface interpolating the surrounding curves, which correspond to the contour-sets on the two planes; see §2. Thus, for the case where no correspondence between contours can be established, this approach offers a possible solution taking also into consideration the shape information encoded in the surrounding curves. Furthermore, enriching each contour with cross-plane tangent distributions, one can handle the global reconstruction problem involving contours on several parallel planes.

The rest of the paper is structured in three sections. Section 2 starts with a coarse description of our method by presenting its key features, namely the quadrilateral tiling of a disc with holes homeomorphic with the sought for surface, the surrounding curve and the Voronoi vertices enclosed therein and, finally, the minimal control cage and its refinement. The section ends with §2.2, providing the algorithm for constructing the initial minimal control cage and the way it is refined for improving the quality of the surface. In §3, we test the performance of the algorithm for two one-to-many problems: a simple radially symmetric 1-3 configuration with circular contours (see Fig. 7) and a complex 1-18 configuration with contours possessing size and shape ranging over a wide spectrum; see Fig. 9. Finally, two additional examples, Figs. 12, 15, the second one involving real medical data, illustrate that the algorithm can be readily used for handling data lying on several planes. The paper ends by summarizing the basic advantages of the proposed approach and discussing future work ( §4). Appendix A outlines the adopted approach for refining the control polygons of the given contours, so that the resulting curves, independently of their parameterization, are within a specified tolerance from the initial contours. This approximation step is required so that the generated T-spline surface approximates, at a user-specified tolerance, the input contours.

Method description

General

As mentioned in the introduction, the paper focuses on developing an automatic algorithm for solving the one-to-many branching problem formulated as below:

• input: a set of planar, simple and closed curves represented as B-Spline curves with individual parameterizations: one base contour C 0,1 on one plane (plane-0), plus multiple contours C 1,µ , µ = 1, ..., m on a parallel plane (plane-1), and a tolerance value limiting the allowed deviation from the input curves; see Fig. 1.

• output: construct a T-spline surface approximating, at tolerance , the input contours on the two planes. The construction of the interpolating/approximating T-spline surface is based on the generation of an appropriate control network, the so-called control cage. For this control cage, we need to establish the minimum number of control vertices, required per contour, and their connectivity so that the topological structure of the surface can be established. The algorithmic procedure that computes these minimum numbers and the corresponding connectivity is described in §2.2.

plane-0 plane-1 c0, 1 c1, 1 c1, 2 c1, 3 c1, 4 c1, 5 c1, 6
The surface we are willing to construct, for the aforementioned one-to-many case, is topologically equivalent to a disc with holes. The external circular boundary of the disc corresponds to the single contour on plane-0, whereas the hole boundaries correspond to the set of contours on the remaining plane-1; see Figs 2 and1. Aiming to construct a branching surface comprising a set of quadrilateral surface patches, we would like to generate a surface decomposition of the disc in Fig. 2 into quadrilateral components according to the following rules:

• If a patch corner touches one of the contours, then one of the patch-edges forming this corner should be a part of the touching contour. As a consequence, patch corners on contours are shared by exactly two neighboring patches.

• The number of corners that are common to three or more than four patches should be limited to the extend possible.

0-1b 0-1a 1-1 c0, 1 c1, 1 c1, 2 c1, 3 c1, 4 c1, 5 A 0-1a c0, 1 c1, 1 c1, 2 c1, 3 c1, 4 c1, 5 B Figure 3: Disc decompositions.
Two disc decompositions that follow the previously mentioned rules are depicted in Fig. 3. In both sketches all bullets correspond to patch corners while curve segments between bullets, both solid and dashed, correspond to patch edges. Specifically, black bullets are corners that lie on contours, gray bullets are corners shared by exactly four patches while the white ones correspond to corners that do not lie on contours and are shared among three or more than four patches. Solid edges have one of their endpoints lying on a contour while dashed ones have no endpoint on a contour and are referred to as bridge-edges. We can also view this decomposition as a graph where bullets correspond to vertices and patch-edges to graph-edges. Using the bridge-edges we can establish an association between different contour parts. For this purpose, we classify bridgeedges in three types as follows: a) simple bridges connecting simple vertices, i.e., vertices of degree 4, b) complex edges connecting complex vertices, i.e. vertices of degree greater than 4 and c) intermediate edges where one vertex is simple and the other complex. Hence, two contour parts belonging to two different contours may be either disassociated, i.e., there exists not a single bridge-edge of any kind lying between them or they are associated via:

1. A complex bridge forming the so-called 1-1 connection, i.e., associated contour parts belong to coplanar contours; see Fig. 3. 2. A simple bridge forming the so-called 0 -1a connection, i.e., contour parts lie on different planes and are associated in a simple way: if we removed the bridge, a single quadrilateral patch could be formed with those two parts opposite to each other; see Fig. 3. 3. An intermediate edge forming the so-called 0 -1b connection, i.e., contour parts lie on different planes but the their association is not simple as the removal of the bridge would not form a quadrilateral patch; see Fig. 3.

The difference between the two alternative decompositions A, B, depicted in Fig. 3, is that in decomposition A all associations of the contour parts of C 1,3 belong to 1 -1 category, while in decomposition B, C 1,3 parts have associations that are of all types. In other words, in the first decomposition C 1,3 has no direct association to the contour on plane-0, which is not the case in the second decomposition; see Fig. 3B.

Our method primarily consists of an automatic process that generates the graph in Fig. 3 and uses it as the initial "control cage" for the interpolating surface. The main components in this process deal with the identification of the contour parts from plane-1 that are associated with parts of the single contour on plane-0 and the identification of the white bullets, i.e., the extraordinary points in the controlling network, which control the remaining inter-contour associations. Furthermore, a preprocessing step is likely to be required for guaranteeing that the generated surface will deviate from the input contours within a user-specified tolerance; see Appendix A.

Appealing to the concept of the surrounding curve we can handle the first issue, i.e., finding the associated parts between the two planes, while the set of vertices of the Voronoi diagram of the contours in plane-1, clipped by the surrounding curve, play the role of the extraordinary points. Specifically, the surrounding curve is based on the convex hull of the contours in plane-1. In its simplest form, introduced in [START_REF] Gabrielides | G 1 -smooth branching surface construction from cross sections[END_REF], it coincides with it and in this case the contour parts from plane-1, that connect to plane-0, are those that are on the convex hull. We have generalized this concept (see §8.2 ibid.) so that, given a user-specified angular tolerance, the surrounding curve is permitted to touch contours that are fully in the interior of the convex hull. The extraordinary points can be used for indicating the saddle points in the generated surface and the Voronoi vertices seem a natural selection for their position. As can be seen from Fig. 3, the exact number of such points is affected from the connections to plane-0 which are dictated by the construction of the surrounding curve as described in §2.2. Once we have established the initial control cage, i.e., identified the number of required control points and their connectivity, the next step involves the replacement of contour points with control points of the contours. Finally, the initial control-cage is enhanced by including all contour control-points and introducing additional edges based on the already established connection-correspondence defined by the bridge-edges. This enriched control cage is used for the final, in our implementation cubic, T-spline surface generation, which is C 2 everywhere with the exception of the neighborhood of extraordinary points, where G 1 continuity is achieved.

Construction of the Control Cage.

In this subsection we present the construction process of the T-spline surfaces control cage. We begin by describing the procedure for generating the minimal control cage, i.e., the cage that utilizes the minimum number of contour-points required per contour in both planes, and their connectivity so that its topological structure materializes the quadrilateral tiling needed for the generation of the branching surface, which is topologically equivalent to a disc with holes, as outlined in the first part of this section; see also Fig. 3 . It is worth noticing that in the below construction we assume that all constructed Voronoi vertices are of degree 3 and every triad of contours corresponds to at most one Voronoi Vertex. The generation of the minimal control cage is accomplished by following the steps described in sequel:

Step 0 -discretization: Create a densely discretized representation of the initial contours belonging on both planes. The discretization is controlled by the value of a discretization parameter.

Step 1 -convex-hull and contour connectors on plane-1: Compute the convex hull (CH) for the contour point-sets on plane-1 and identify the line segments, henceforth referred to as CH ContourConnectors, that belong to the convex hull and connect consecutive contour point-sets; see dashed line segments in Fig. Step -surrounding curve & correspondence to plane-0 contour: Assemble the Surround-ingCurve (SC) as the polygonal line that comprises the ContourConnectors, the segments between the separated vertices of Step 4 and the segments joining the endpoints of the Con-tourConnectors incident at the same contour point-set; see Fig. 5a. Using the convex-hull projection, see [START_REF] Ekoule | A triangulation algorithm from arbitrary shaped multiple planar contours[END_REF], generate a polar parameterization of the contour point-set on plane-0 and the SC on plane-1. Next, create line segments connecting the end-points of the ContourConnectors with the points of the contour of plane-0 that share the same polar parametrization; see Fig. 5b. Furthermore, identify the points on the contour point-set on the plane-0, that correspond to the mid-points of the ContourConnectors on plane-1; see dashed lines in Fig. 5b.

Step -linking Voronoi vertices to contours on plane-0 and plane-1: Perform a vertical translation of each Voronoi vertex using the input parameters v immersion i, thus positioning each Voronoi vertex between the two levels. Next, for each ContourConnector, create a line segment connecting the corresponding Voronoi vertex, i.e., the one whose circle touches both connected contours, with the points on the contour of plane-0, identified in Step 5, which correspond to the ContourConnector mid-points; see Fig. 5c.

Step -bridge-edges: Create bridge-edges as follows:

1 -1: Generate line segments connecting neighboring Voronoi vertices, defined as those Voronoi vertices that are associated to the same two contours; see thick green lines in Fig. 5d. 0 -1b: For each ContourConnector, create two line segments, connecting the corresponding Voronoi vertex with the segments connecting plane-0 with plane-1 and emanating from the end-points of the ContourConnector, at positions defined by the parameter v immersion i ; see thick blue lines in Fig. 5d. 0 -1a: Create the line segments connecting the end-points of the 0 -1b edges that lie on lines emanating from the same contour of plane-1; see thick red lines in Fig. 5d. Step 8 -minimal control cage: Connect by linear segments the identified points on each contour (in both planes); see thick solid lines in Fig. 5e. These segments along with the Voronoi vertices and the complete connectivity of points and Voronoi vertices described in the previous steps constitute the minimal control cage; see Figs 5e,5f.

Having created the minimal control cage, we proceed with refining the cage in order to guarantee that the final T-spline surface deviates from the given contours within a tolerance , regardless of the parameterizations employed in its construction. To this end, we first employ the algorithm described in Appendix A for a given tolerance and produce control polygons for all given contours; see thick dashed lines in Fig. 6a. If the generated control points of a contour are fewer than implied by the minimal control cage we perform adequate knot-insertions. We then relocate the contour-points of the minimal control cage to their nearest control points of the individual contours, leading to a modified control cage and an explicit correspondence map between the control points of different contours via the bridge-edge concept; see Fig 6a . This correspondence map is used for refining the cage, since the parts requiring refinement can be easily identified. The refining process comprises the construction of the polygonal lines connecting the remaining control points, the modification of the bridge edges in order to encapsulate the shape of the connected polygonal lines and finally, the creation of linear T-edges between corresponding polygonal lines and bridge edges; see red, green and blue thick lines in Fig. 6a corresponding to 1 -1, 1 -0a and 1 -0b associations respectively. The final control cage, depicted in Fig. 6b, is used to construct the T-spline branching surface depicted in Fig. 6c.

Examples

Our method, described in §2, has been implemented as a plug-in to Rhinoceros [START_REF] Mcneel | Rhinoceros, design, model, present, analyze, realize[END_REF] modeling environment with the aid of Autodesk 's T-Splines Plug-In [START_REF]Autodesk Inc, Autodesk T-Splines Plug-in for Rhino[END_REF] for Rhinoceros . Our plug-in is coded in C# using the RhinoCommon SDK 1 . All examples presented in this work have been produced using the aforementioned implementation.

Our first example involves a simple data-set with one contour on plane-0 and three contours on plane-1. All contours are circles and the data-set is radially symmetric as can be seen from Figures 7 and8. The generated minimal (initial) and final control cages are depicted in Fig. 7a and 7b, while Fig. 7d presents a renders view of the final T-Spline surface based on the final control-cage shown in Fig. 7c. Finally, in Fig. 8, the Gaussian curvature for the T-spline surfaces based on the minimal (initial) and final control cages are depicted. For both cases it is easy to observe that the radial symmetry of the data-set is preserved in our construction. In this example, we have a single extraordinary point in the construction and the overall number of vertices for the minimal control cage is 43 while this number rises to 140 when we require that the circular contours are approximated with a tolerance of 1%. The conversion2 of the latter T-spline surface to a NURBS representation yields a multi-patch surface with 1416 control points.

In our second example, we are demonstrating the output of our approach for the case of a more complex data-set with 18 contours on plane-1 of varying size and shape; see Fig. Figure 10 depicts the generated T-spline surface using the initial / minimal control cage. Its enrichment with the total number of control points used in the representation of contours leads to the final control cage which generates the T-spline surface depicted in Fig. 11. This final surface guarantees that distance of the boundary edges of the T-spline surface from the corresponding input contours are below a specified tolerance, which, in this case, is below 0.05% of contour's length. The number of control points in this final representation is 741 and the conversion to a NURBS yields 126 surface patches with an overall number of control points rising to 28942. Finally, the number of extraordinary points for this example is 24.

In our final examples, we will demonstrate the usage of the developed approach for the generation of T-spline surfaces interpolating / approximating contours in multiple planes. The first of these examples involves 3 planes as can be seen in Fig. 12. The overall number of contours is 8 with 1 contour lying at each boundary plane (plane-1 and plane-1*) and the remaining 6 contours at the intermediate plane-0. We apply our construction process twice, generating one T-spline surface, T 1 , between the contours on planes 1 and 0 and one surface, T 1 * , between planes 0 and 1*. The control cages for those two surfaces are depicted in Fig. 14a. Employing the same approximating tolerance, the boundary edges of those T-spline surfaces on plane-0 will coincide. As a second step, in our construction for this particular example, we modify the second row of control points in both surfaces so that, for each boundary control point b i on plane-0, the vectors b i,1 -b i and b i,1 * -b i are collinear, where b i,1 , b i,1 * denote the previous and next to b i control points that belong to T 1 and T 1 * , respectively; see Fig. 14b. Employing this modification for all the boundary control points on plane-0, first-order geometric continuity is secured between the two surfaces, as can been seen in Fig. 13. The extraordinary points in this final construction are 8 overall; 4 in each T-spline surface.

As our last example, we use the 44 contours shown in Fig. 15a coming from a mesh model of the eighth tooth retrieved from the Aim@Shape shape repository [START_REF]aim@shape, Aim@Shape Repository[END_REF]. The original model is a surface mesh in VRML format, which is depicted in Fig. 15d. The contours were acquired by initially intersecting the initial mesh with parallel planes and subsequently approximating them with cubic B-Splines. The one-to-many cases exhibited in this example are three, denoted with think solid lines in Fig. 15a. Their resulting T-spline control cages are shown in Fig. 15b, while the reconstructed surface for the whole tooth is depicted in Fig. 15c. Obviously, the surface approximation quality is dependent on the selection and approximation of the contours set, however, with this rather rough set we obtained a maximum deviation, between the initial and reconstructed surface, of less than 1% of tooth's width. The required number of control points for all three branching areas is around 1000, while more than 8000 would be required for acquiring the same approximation level using a multi-patch NURBS surface approach.

Conclusions & Future work

In this work, we have presented an automatic construction process for the generation of a smooth T-spline surface that interpolates / approximates simple planar curves (contours) that lie on parallel planes being C 2 everywhere except from the neighborhood of a finite set of points, where its smoothness degrades to G 1 level. Specifically, we have demonstrated the automatic generation of the T-spline control cage (T-mesh) that defines the T-spline surface for the special case of the "one-to-many" problem between two planes, where multiple contours lie on one plane and only one on the other. In our construction we have assumed that Voronoi vertices are of degree 3 and that for every triplet of contours there's a maximum of one Voronoi vertex. These assumptions are used for simplifying the construction of complex and intermediate bridge-edges. However, if we incorporate the actual Voronoi edges in our construction, the requirement for these limiting assumptions can be eliminated.

The proposed construction can be also utilized for the case of the more general "many-to-many" problem, where multiple contours lie on both planes. This can be accomplished by the introduction of a single contour on an intermediate plane. The introduced contour is not an arbitrary contour but is generated by intersecting the surface interpolating the surrounding curves from the two planes with the intermediate plane. Hence, the introduced contour captures the shape information from the two planes and allows us to split the "many-to-many" problem in two "one-to-many" problems. Obviously, this solution is limited by the assumption that all contours converge to a single one lying on the introduced plane, which however is a reasonable assumption when no contours' correspondence can be established. A more elaborate procedure must be followed if we want to avoid this convergence to a single contour which could be a topic for future research. Figure 14: The control cage of the T-spline surface interpolating contours on 0 and 1 planes is black, while the one corresponding to planes 0 and 1* is gray. Tangency handles are modified locally about plane-0 so that we can guarantee tangential continuity.

Furthermore, we have demonstrated the use our construction algorithm for the case of 3 parallel planes (see example 3 in §3) where the introduction of cross-plane tangent distributions on along the contours of the intermediate planes permits the extension of our approach to more than 2 planes. Finally, our approach is applied for reconstructing a tooth surface based on a rather small number of 44 contours. The resulting reconstructed surface is considerably close to the initial tooth-mesh as already discussed in §3.

Currently, due to certain simplifying assumptions, the implementation of our algorithm is restricted to contours that lie on parallel planes. This limitation could be overcome by the introduction of a different way for separating contour sets and a more general way for performing the Voronoi vertices' immersion in step 6 of our construction process. Furthermore, in our current implementation we have not taken any precautions for guaranteeing that the generated T-spline surface is analysis suitable, however, if required, we may easily accomplish this by appropriate refinement in the neighborhood of extraordinary points. Some of the simplifying assumptions in our current implementation are due to the employed software package, i.e., Autodesk 's T-Splines Plug-In, and the restrictions set by its use. Hence, it would be beneficial for the enhancement of the method and its implementation, if we could access the API of the specific software package or alternatively use a T-spline software library that allows access to internals of the parameterization.

Finally, it is among our aims to pursue the further development of this construction method for trivariate T-splines that would enable volumetric meshing and open a new window for possible applications. A possible way towards this aim may come out of the combination of the T-spline surface representation and the medial axis transform.

Appendix A. Controlling the deviation of the T-spline surface from the contour curves

Let us assume that all of our input contours are represented by closed cubic B-Spline curves and the generated surface is a bicubic T-spline surface with control polygons of boundary edges that coincide with those of the input contours. In this case, if we could impose an appropriate parameterization for each boundary T-spline surface edge, the generated surface would interpolate the input contours. Unfortunately, imposing the appropriate parameterization for all edges would explode the complexity of both the control-cage generating process and the resulting T-spline surface. In this work, we employ a simple uniform parameterization defined by the T-spline construction process; see §2.2. Hence, although the control polygons of the contours coincide with the control polygons of the T-spline boundary edges, the parameterizations of the boundary edges are generally not the ones used in the input contours and as a result, the T-spline surface does not interpolate them. Two B-spline curves with different parameterizations, which however share the same degree and control polygon are confined to lie within the same area defined by the union of the local convex-hulls generated by corresponding control points; see Fig. Based on this observation we can devise a simple process which can generate a control polygon that when shared by two curves, with arbitrary parameterizations, will confine them to be within a given distance from each other. Let's assume that two cubic B-Spline curves, c 1 (t 1 ), c 2 (t 2 ), t 1 ∈ [a, b] and t 2 ∈ [c, d] with different knot vectors use the same control polygon. The distance between them, measured by the Fréchet metric, is defined to be: respectively. Let us also assume that the union U of their local convex-hulls is defined by the polygonal curves B i and B e . As both curves lie in U , it is obvious that δ F (c 1 , c 2 ) ≤ δ F (B i , B e ). Furthermore, for every possible reparameterizations t i and t j of c 1 and c 2 the above inequality still holds. Thus, if δ F (B i , B e ) ≤ the distance of any pair of cubic B-Spline curves with the same control polygon but arbitrary parameterizations will be less than , i.e., δ F (c 1 (t i ), c 2 (t j )) ≤ .

δ F (c 1 (t 1
As the computation of the exact Fréchet distance is a fairly complicated procedure, we use the discrete Fréchet distance (δ dF ) and the algorithmic procedure for its computation suggested by Eiter and Mannila [START_REF] Eiter | Computing discrete Fréchet distance[END_REF]. Since, for all polygonal curves P and Q: δ F (P, Q) ≤ δ dF (P, Q), we may write that: δ F (c 1 (t i ), c 2 (t j )) ≤ δ F (B i , B e ) ≤ δ dF (B i , B e ).

Obviously, if δ dF (B i , B e ) ≤ ⇒ δ F (c 1 (t i ), c 2 (t j )) ≤ . Hence, the procedure for assuring that the input contour c j and the corresponding T-spline boundary-edge will lie within an distance has as follows:

Figure 1 :

 1 Figure 1: Input data for the one-to-many problem.

Figure 2 :

 2 Figure2: Hollowed disc corresponding to interpolating surface. Surface borders correspond to contours on plane-0 (C 0,1 ) and plane-1 (C 1,i , i = 1, . . . , 5).

4a. Step 2 -

 2 Voronoi triangles identification: Compute the Delaunay triangulation of the pointsets on plane-1 and identify the Voronoi triangles, i.e., the triangles connecting points from 3 distinct contour point-sets; see solid thick triangles in Fig.4b. Step 3 -modification of contour connectors: Replace each CH ContourConnector with the internal sides of the corresponding Voronoi triangle, if the triangle's internal angle ω satisfies ω ≥ α, where α is an input angular parameter, and generate the set of ContourConnectors; see Figs. 4b-c. This procedure is carried out iteratively until all Voronoi triangles have an internal angle smaller than α; see Figs. 4b-d. The triangles that have been used in this procedure are removed from the construction process. Step 4 -connecting Voronoi vertices to contours: For each remaining Voronoi triangle generate the corresponding Voronoi vertex. Moreover, for each pair of contours that is connected with two ContourConnectors, introduce an additional vertex on the corresponding Voronoi edge, which will act as an additional Voronoi vertex. Next, for each Voronoi vertex, generate a connecting line segment between the vertex and the corresponding point on each contour; see Fig. 4e. Finally, if two ContourConnectors have a common end-vertex, this vertex is replaced by two new ones on the corresponding contour using an appropriate separation parameter and the final set of ContourConnectors is produced; see dashed lines in Fig. 4f.

Figure 4 :

 4 Figure 4: Minimal control cage construction (steps 0-4): working on plane-1.

Figure 5 :

 5 Figure 5: Minimal control cage construction (steps 5-8).

  Figure 6

  9. As can be easily seen from this figure, 4 contours (C 1,5 , C 1,8 , C 1,9 and C 1,17 ) are non-convex with C 1,8 exhibiting a very abrupt shape change.

Figure 7 :

 7 Figure 7: Control cages of a radial-symmetric data-set with 4 input contours: 1 on plane-0 and 3 on plane-1.

Figure 8 :

 8 Figure 8: Gaussian curvature graph for the initial (left) and final (right) surface using cages shown in 7b and 7c, respectively.

Figure 9 :

 9 Figure 9: Data-set with 19 input contours: 1 on plane-0 and 18 on plane-1.

Figure 10 :

 10 Figure 10: Rendered view of the initial T-spline surface along with the input contour-set.

Figure 11 :

 11 Figure 11: Rendered view of the final T-spline surface along with the input contour-set.

Figure 12 :

 12 Figure 12: Data-set with contours on 3 parallel planes: plane-1, plane-0 and plane-1* Figure 13: Rendered view of the final T-spline surface. Contours are depicted in dark-gray color.

Figure 15 :

 15 Figure 15: Model of tooth No.8, retrieved from the Aim@Shape shape repository

Figure A. 16 :

 16 Figure A.16: Two cubic B-Spline curves lying in the union of their local convex hulls: (b 1 b 4 b 2 ) ∪ (b 4 b 5 b 2 b 3 ) ∪ (b 4 b 6 b 3 ) ∪ (b 4 b 7 b 6 b 5 ).Dashed curve's knot vector is (0, 0.872, 0.962, 0.983, 1) while the solid one uses: (0, 0.046, 0.656, 0.702, 1).

Figure A. 17 :

 17 Figure A.17: (a) A cubic B-Spline along with its initial control polygon and local-convex-hull union U 1 , colored gray. (b) The same B-Spline curve after uniform knot insertion along with the corresponding local-convex-hull union U 2 . (c) After an additional uniform refinement, ∂U 3 is defined by two polygonal boundary curves.

), c 2 (t 2 d (c 1

 221 (α(t)) , c 2 (β(t))) ,where d is a 2D point metric (e.g., Euclidean distance) and t 1 = α(t), t 2 = β(t) are continuous, non-decreasing functions mapping [0, 1] to [a, b] and [c, d],

1 .

 1 Set k = 0 2. Compute the union of the local convex hulls (U k ) for the control polygon of c jk . 3. If U is not defined by exactly two separate polygonal lines, perform uniform knot insertion to the contour until U has the required property. 4. Assuming U is defined by polygonal curves B ik and B ek , compute δ dF (B ik , B ek ). 5. While δ dF (B ik , B ek ) ≥ • identify the area of c jk where δ dF (B ik , B ke ) attains its value. Perform knot insertion in the identified segment of c jk and get the new c i,k+1 contour • Set k=k+1 • Compute the new U k and the corresponding B ik , B ek 6. Return the final c jk .

RhinoCommon is a cross-platform .NET plug-in Software Development Kit.

The conversion is performed using the corresponding functionality of Autodesk 's T-Splines Plug-In.