Ioannis Z Emiris

Konstantinos Gavriil

Christos Konaxis

Interpolation of syzygies for implicit matrix representations

Keywords: matrix representation, syzygies, implicitization, point cloud, interpolation, space curve, triangular surfaces

We examine matrix representations of curves and surfaces based on syzygies and constructed by interpolation through points. They are implicit representations of objects given as point clouds. The corresponding theory, including moving lines, curves and surfaces, has been developed for parametric models. Our contribution is to show how to compute the required syzygies by interpolation, when the geometric object is given by a point cloud whose sampling satisfies mild assumptions. We focus on planar and space curves, where the theory of syzygies allows us to design an exact algorithm yielding the optimal implicit expression. The method extends readily to surfaces without base points defined over triangular patches. Our Maple implementation has served to produce the examples in this paper and is available upon demand by the authors.

Introduction

In manipulating curved geometric objects, it is essential to possess robust algorithms for changing representation. This paper considers the representation of an object by point samples and offers an algorithm which constructs a matrix representation of implicit curves and surfaces by interpolating the algebraic syzygies of the (unknown) parametric equations through points. Point samples are very common in industrial applications; this paper proposes a model for using them in the framework of syzygies.

A matrix representation of an implicit object is a single matrix, generically of full rank, which represents the object in the sense that its rank drops precisely when evaluated at a point lying on the object. Matrix representations are quite robust, since they do not require computation of the implicit equation; instead, they reduce geometric operations on the object to linear algebra. In general, existing approaches to implicitization include Gröbner bases, resultants, moving lines/curves and surfaces, µ-bases and approximation complexes, as well as a number of interpolation techniques. Today, moving lines/curves and surfaces, and µ-bases seem to offer very competitive methods since they provide the veracity of algebraic approaches without the high complexity of Gröbner bases nor the problems due to base points when using resultants. Moving curves and surfaces have been used to construct matrix representations of implicit objects, and this is the premise of our work.

The theory of syzygies, including moving lines, curves and surfaces, has been developed for parametric models; it is sketched in the subsequent sections. Our contribution is to show how to compute the required syzygies by interpolation, when the input curve or surface is given by a point cloud whose sampling satisfies mild assumptions. No information on the parametric representation of the object is given, but the parametric expressions could be obtained from the algorithm's output. However, our goal is a robust implicit matrix representation, and we focus on matrices constructed only by linear syzygies. We illustrate our algorithms for planar and space curves as well as triangular surfaces, all without base points. Our Maple implementation is available upon demand by the authors.

Let us describe the input in the case of curves in an ambient space of arbitrary dimension n ≥ 2; it shall be generalized in the sequel to surfaces. We assume that the curve admits some (unknown) affine rational parameterization φ : R → R n ; planar and space curves correspond to n = 2 and n = 3, respectively. The input is a parametric set of points. This pointset is defined as a sequence of vectors (τ k ; X k ) such that φ(τ k ) = X k , for all k = 1, 2, . . . , where τ k ∈ R, and X k ∈ R n . In particular, X k ∈ R 2 or X k ∈ R 3 , depending on whether we study planar or space curves. For sampling such a set of points, φ may be an arc-length parameterization and the vectors can be obtained by a scanner capable of measuring the distance it has covered when moving on the curve. In practice, this assumption can be satisfied when the scanner is equipped with a GPS system. When given a set of points that is dense enough, the distances between consecutive points can be used to approximate the arc length of the curve.

A related model for point clouds is considered in [START_REF] Floater | Parameterization for curve interpolation[END_REF]. This paper is organized as follows: Section 2 overviews previous work, whereas Section 3 contains some background in the theory of syzygies, and develops general tools required in the sequel. Section 4 describes our method for interpolating syzygies when the input is given as a set of parametric points defining a planar or space curve. In Section 5 we extend the method to the case of triangular surfaces, given as a parametric pointset. We conclude with future work and open questions.

Previous work

This section discusses the main existing approaches to implicitization, with an emphasis on methods constructing matrix representations of implicit objects. Besides these methods, Gröbner bases offer a powerful and complete approach but suffer from high complexity and numerical instability.

Resultants, and their matrix formulae, have been used to express the implicit surface equation, e.g., in [START_REF] Manocha | The implicit representation of rational parametric surfaces[END_REF], under the assumption of no base points. The most direct method to reduce implicitization to linear algebra is to construct a square matrix M , indexed by all possible monomials in the implicit equation (columns) and different values (rows) at which all monomials get evaluated. Then the vector of coefficients of the implicit equation is in the kernel of M . This idea has been extensively used, e.g. in [Dok01, EKKL13, EKK15, SY08]. The method, as introduced in [EKKL13, EKK15], exploits sparse resultant theory so as to predict the monomials in the implicit equation and thus build the interpolation matrix. It handles objects with base points.

A modern method for representing implicit equations by matrices was introduced by Sederberg and his coauthors when they rediscovered the theory of syzygies in the context of computer science [START_REF] Sederberg | Curve implicitization using moving lines[END_REF][START_REF] Sederberg | Implicitization using moving curves and surfaces[END_REF][START_REF] Sederberg | Implicitizing rational curves by the method of moving algebraic curves[END_REF]. Let us take the example of planar curves without base points, parameterized by the homogeneous polynomials (f 1 (s : t) : f 2 (s : t) : f 3 (s : t)), all of same degree d. The main idea is to define a moving line in P 2 as

h 1 (s : t)x + h 2 (s : t)y + h 3 (s : t)z = 0, (1)
where x, y, z are homogeneous coordinates in P 2 and h i (s : t) ∈ C[s, t], i = 1, 2, 3, are homogeneous polynomials of same degree. The moving line follows the curve if

3 i=1 h i (s : t)f i (s : t) = 0, for all (s : t) ∈ P 1 . (2)
Algebraically, the triplet (h 1 , h 2 , h 3 ) of homogeneous polynomials h i , or, equivalently, the moving line (1), is a (linear) syzygy on the polynomials f i . It is known, see e.g. [START_REF] Sederberg | Curve implicitization using moving lines[END_REF][START_REF] Cox | Equations of parametric curves and surfaces via syzygies[END_REF], that there are d independent moving lines of degree d -1 that follow the curve. Using these moving lines it is possible to construct a d × d matrix whose determinant is a multiple of the implicit equation, see Proposition 1.

In the next section we provide a comprehensive discussion on syzygies. For now, let us recall that in the case of surfaces without base points, one may also construct a square matrix whose determinant is a power of the implicit polynomial [START_REF] Cox | On the validity of implicitization by moving quadrics for rational surfaces with no base points[END_REF], by using d moving planes and (d 2d)/2 moving quadrics, all of degree d -1, see Subsection 3.4.

If we allow orthogonal matrices, it suffices to work with linear syzygies, and this is the main approach adopted in this work. In general, one defines the notion of critical degree ν 0 , see Proposition 5, which corresponds to the degree of the linear syzygies required to define an orthogonal matrix M ν (φ) that satisfies the following property [START_REF] Busé | Matrix-based implicit representations of rational algebraic curves and applications[END_REF]: for any point p ∈ P 2 in the case of planar curves or, respectively, p ∈ P 3 in the case of space curves or surfaces, the rank of M ν (φ) evaluated at p drops if and only if p belongs to the algebraic closure of Im(φ). The critical degree is, in general, at least as large as the regularity of the map sending tuples of polynomials to combinations generalizing those in expression (2). In particular, the critical degree in the case of planar and space curves without base points is d -1, and for triangular surfaces it is 2(d -1).

The matrices indirectly represent implicit objects and allow for geometric operations, such as surface-surface intersection [START_REF] Busé | The surface/surface intersection problem by means of matrix based representations[END_REF] and, more recently, ray shooting [START_REF] Shen | A line/trimmed NURBS surface intersection algorithm using matrix representations[END_REF], to be reduced to linear algebra. Their advantage is that the matrices are much smaller than interpolation matrices, and allow for inversion by an eigenproblem on these matrices. They also simplify in the presence of base points while other methods become more complicated. On the other hand, their construction is a two-step process of matrix operations. Moreover, they are symbolic with entries linear polynomials in the implicit variables.

Basic tools

This section uses known results in the theory of syzygies to develop certain tools needed for stating our algorithms in subsequent sections. In particular we shall relate the degree of a given grading of the syzygy module to its dimension. For a comprehensive survey on the subject, we refer the interested reader to [START_REF] Cox | Equations of parametric curves and surfaces via syzygies[END_REF][START_REF] Cox | Curves, surfaces, and syzygies[END_REF].

Planar curves

Consider the (homogeneous) parameterization of a planar curve C:

φ : P 1 → P 2 : (s : t) → (f 1 (s : t) : f 2 (s : t) : f 3 (s : t)) , (3)
where f i ∈ C[s, t] are homogeneous of the same degree d, and assume that gcd(f 1 , f 2 , f 3 ) = 1, i.e. φ has no base points.

Consider a syzygy (h 1 , h 2 , h 3 ), where h i (s : t) ∈ C[s, t], i = 1, 2, 3, are homogeneous polynomials of same degree, as in (2):

3 i=1 h i (s : t)f i (s : t) = 0, for all (s : t) ∈ P 1 ,
This is a linear syzygy on the polynomials f i . The common degree of the h 1 , h 2 , h 3 is known as the degree of this syzygy. The set of all syzygies is denoted by Syz(f 1 , f 2 , f 3 ), and has the structure of a graded module. By fixing a degree ν ≥ 0, we can consider the set of syzygies of degree ν, denoted by Syz(f 1 , f 2 , f 3 ) ν , which is known to be a finite dimensional C-vector space. One can compute a basis L 1 , . . . , L Nν of this vector space by solving a linear system, where N ν denotes the basis cardinality.

We identify each L j = (h

(j) 1 , h (j) 2 , h (j)
3 ) with its moving line and we develop it in terms of the s, t as follows:

L j := 3 k=1 h (j) k x k = ν i=0 Λ i,j (x, y, z)s i t ν-i , j = 1, . . . , N ν , (4)
where Λ i,j (x, y, z) is a linear polynomial in C[x, y, z]. Let M ν (φ) be the (ν + 1) × N ν matrix, whose jth column contains the coefficients Λ i,j (x, y, z) of L j in (4).

A fundamental result here is the following, showing that there are d independent moving lines of degree d -1 that follow φ. The entire module of syzygies Syz(f 1 , f 2 , f 3 ) is a free module of rank 2. Let P = (P 1 , P 2 , P 3 ),

Q = (Q 1 , Q 2 , Q 3 ) be its generators of degrees µ 1 ≤ µ 2 , respectively. It is known that µ 1 + µ 2 = d.
The P, Q are called a µ-basis of Syz(f 1 , f 2 , f 3 ). Thus, we write any syzygy in Syz(f 1 , f 2 , f 3 ) ν as a polynomial combination, for homogeneous p, q ∈ C[s, t], namely:

pP + qQ, where deg (p) = ν -µ 1 and deg (q) = ν -µ 2 , . (5)
If we identify P, Q with their moving lines, i.e., P = P 1 x + P 2 y

+ P 3 z, Q = Q 1 x + Q 2 y + Q 3 z,
then the Sylvester resultant of P, Q gives the implicit equation F of C:

Res(P, Q) = c • F deg(φ) ,
where c, F, deg(φ) are as in Proposition 1.

We now employ (5) to compute a basis of Syz(f 1 , f 2 , f 3 ) ν and its dimension, as ν varies. Recall d is the homogeneous degree of f 1 , f 2 , f 3 . The following lemma essentially appears in [CSC98, Cor.2,p. 811]. For the convenience of the reader we include a self-contained and simple proof.

Lemma 2. We distinguish the following cases for the degree

ν of Syz(f 1 , f 2 , f 3 ) ν : (a) ν ≤ µ 1 -1. Then dim Syz(f 1 , f 2 , f 3 ) ν = 0. (b) µ 1 -1 ≤ ν ≤ µ 2 -1. Then dim Syz(f 1 , f 2 , f 3 ) ν = ν -µ 1 + 1. (c) ν ≥ µ 2 -1. Then dim Syz(f 1 , f 2 , f 3 ) ν = 2ν -d + 2.
Note that the intervals of the three cases share their endpoints, hence the overall piecewise linear curve is continuous. The lemma generalizes the fundamental result that dim

Syz(f 1 , f 2 , f 3 ) d-1 = d, from Proposition 1.
Proof. From equation ( 5) we get the basis of Syz(f 1 , f 2 , f 3 ) ν , for general ν:

B = {s i t ν-µ 1 -i P | 0 ≤ i ≤ ν -µ 1 } ∪ {s i t ν-µ 2 -i Q | 0 ≤ i ≤ ν -µ 2 }. (6)
Then, the lemma follows straightforwardly by computing the cardinality of B for each case (a)-(c). In particular we have:

(a) If ν ≤ µ 1 -1, then B = ∅ because any non-trivial polynomial combination of P, Q has total degree ≥ µ 1 . Hence dim Syz(f 1 , f 2 , f 3 ) ν = 0. (b) If µ 1 ≤ ν ≤ µ 2 -1, then B = {s i t ν-µ 1 -i P | 0 ≤ i ≤ ν-µ 1 } and |B| = dim Syz(f 1 , f 2 , f 3 ) ν = ν -µ 1 + 1. If ν = µ -1 this formula yields correctly 0. (c) If ν ≥ µ 2 , then B is as in (6), containing both multiples of P and Q, hence |B| = dim Syz(f 1 , f 2 , f 3 ) ν = 2ν -d + 2, since d = µ 1 + µ 2 . At ν = µ 2 -1, the formula yields µ 2 -µ 1
, which is also obtained at this point by the formula of case (b).

The lemma is summarized in Figure 1.

Space curves

Consider the space curve parameterized homogeneously as

φ : P 1 → P 3 : (s : t) → (f 1 (s : t) : f 2 (s : t) : f 3 (s : t) : f 4 (s : t)), (7)
where d is again defined as the homogeneous degree of the polynomials f i (s, t), i = 1, . . . , 4. Suppose gcd(f 1 , f 2 , f 3 , f 4 ) = 1, i.e. there are no base points.

The module of syzygies Syz(f 1 , f 2 , f 3 , f 4 ) is a free module of rank 3. Let P = (P 1 , P 2 , P 3 , P 4 ),

Q = (Q 1 , Q 2 , Q 3 , Q 4 ), R = (R 1 , R 2 , R 3 , R 4
) be its generators and µ 1 ≤ µ 2 ≤ µ 3 be their degrees respectively. It is known that

µ 1 + µ 2 + µ 3 = d.
The P, Q, R are called a µ-basis of Syz(f 1 , f 2 , f 3 , f 4 ). We can write any syzygy in Syz(f 1 , f 2 , f 3 , f 4 ) ν as a polynomial combination for homogeneous polynomials p, q, r ∈ C[s, t]:

pP + qQ + rR, where deg (p) = ν -µ 1 , deg (q) = ν -µ 2 , deg (r) = ν -µ 3 . (8)
Identifying P, Q, R with their moving lines, i.e., P = P 1 x + P 2 y + P 3 z

+ P 4 w, Q = Q 1 x + Q 2 y + Q 3 z + Q 4 w, R = R 1 x + R 2 y + R 3 z + R 4 w
, and forming the Sylvester resultant of every pair of P, Q, R gives one implicit equation of a surface containing curve C; the latter is thus defined set-theoretically as the intersection of 3 surfaces.

We can now relate the dimension of

Syz(f 1 , f 2 , f 3 , f 4 ) ν to ν. µ 1 -1 µ 1 µ 2 -1 µ 2 d -1 d µ 2 -µ 1 + 2 µ 2 -µ 1 ν N ν ν - µ 1 + 1 2 ν -d + 2 d ν + 1
Figure 1: The graph of the dimension N ν of Syz(f 1 , f 2 , f 3 ) ν with respect to ν. The dashed red line intersects the graph at point (d -1, d) corresponding to the critical degree.

Lemma 3. We distinguish the following cases for the degree

ν of Syz(f 1 , f 2 , f 3 , f 4 ) ν : (a) ν ≤ µ 1 -1. Then dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 0. (b) µ 1 -1 ≤ ν ≤ µ 2 -1. Then dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = ν -µ 1 + 1. (c) µ 2 -1 ≤ ν ≤ µ 3 -1. Then dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 2ν -µ 1 -µ 2 + 2. (d) µ 3 -1 ≤ ν. Then dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 3ν -d + 3.
The intervals of subsequent cases share their endpoints, hence the overall piecewise linear curve is continuous.

Proof. From equation (8) we get the basis of Syz(f 1 , f 2 , f 3 , f 4 ) ν for general ν:

B = {s i t ν-µ 1 -i P | 0 ≤ i ≤ ν-µ 1 }∪{s i t ν-µ 2 -i Q | 0 ≤ i ≤ ν-µ 2 }∪{s i t ν-µ 3 -i W | 0 ≤ i ≤ ν-µ 3 }.
(9) Then, the lemma follows straightforwardly by computing the cardinality of B for each case (a)-(d). In particular we have:

(a) If ν ≤ µ 1 -1, then B = ∅ because any non-trivial polynomial combination of P, Q, W has total degree ≥ µ 1 . Hence dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 0. (b) If µ 1 ≤ ν ≤ µ 2 -1, then B = {s i t ν-µ 1 -i P | 0 ≤ i ≤ ν-µ 1 } and |B| = dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = ν -µ 1 + 1. If ν = µ 1 -1 the formula yields correctly 0. µ 1 -1 µ 2 -1 µ 2 -µ 1 ν N ν ν -µ1 + 1 2 ν - µ 1 - µ 2 + 2 3 ν - d + 3 µ 3 -1 2µ 3 -µ 2 -µ 1 2 ( ν + 1 ) d -1 2d Figure 2: The graph of the dimension N ν of Syz(f 1 , f 2 , f 3 , f 4 ) ν with respect to ν. The dashed red line intersects the graph at point (d -1, 2d) corresponding to critical degree ν 0 = d -1. (c) If µ 2 ≤ ν ≤ µ 3 -1, then B = {s i t ν-µ 1 -i P | 0 ≤ i ≤ ν -µ 1 } ∪ {s i t ν-µ 2 -i Q | 0 ≤ i ≤ ν -µ 2 } and |B| = dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 2ν -µ 1 -µ 2 + 2. If ν = µ 2 -1 the formula yields µ 2 -µ 1
, which is also obtained from the formula of case (b).

(d) If ν ≥ µ 3 , then B is as in (9) and |B| = dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 3ν -d + 3. If ν = µ 3 -1,
the formula yields 2µ 3µ 1µ 2 , which agrees with the value of the formula in case (c) at this point.

Figure 2 summarizes the lemma.

General Curves

We may unify and generalize the previous discussion by considering curves in P n , for any ambient dimension n ≥ 2, parameterized homogeneously as

φ : P 1 → P n : (s : t) → (f 1 (s : t) : . . . : f n (s : t)), ( 10
)
where d is the homogeneous degree of the polynomials f i (s, t), i = 1, . . . , n. By the Hilbert Syzygy Theorem, the syzygy module Syz(f 1 , . . . , f n ) is free of rank n, and in particular it has a µ-basis [CSC98, Thm.1]. Let µ 1 ≤ • • • ≤ µ n , be the degrees of the polynomials in the µ-basis of the module. Hence, the previous discussion extends to this case as well.

The derived formulae for N ν in the two Lemmas above can be unified and generalized into a piecewise linear formula with n nontrivial pieces, where the equation of the k-th segment is

N ν = k i=1 (ν -µ i + 1), for µ k -1 ≤ ν ≤ µ k+1 -1, k = 1, 2, . . . , n -1, or µ n -1 ≤ ν.
Of course N ν = 0 for ν ≤ µ 1 -1.

Triangular surfaces

The theory of moving lines generalizes to surfaces in P 3 . Let us focus on the case of surfaces without base points, parameterized by

φ : P 2 → P 3 : (s : t : u) → (f 1 (s : t : u) : f 2 (s : t : u) : f 3 (s : t : u) : f 4 (s : t : u)) , (11)
where all f i are homogeneous of degree d. These are known as triangular surfaces.

The analogue of a moving line of degree ν in P 3 is a moving plane:

h 1 (s : t : u)x + h 2 (s : t : u)y + h 3 (s : t : u)z + h 4 (s : t : u)w = 0, ( 12
)
where deg(h i ) = ν, i = 1, . . . , 4.

A moving quadric of degree ν is defined as:

h 1 x 2 + h 2 y 2 + h 3 z 2 + h 4 w 2 + h 5 xy + h 6 xz + h 7 xw + h 8 yz + h 9 yw + h 10 zw = 0, (13)
where where deg(h i (s : t : u)) = ν, i = 1, . . . , 7.

There are d linearly independent moving planes of degree d -1 that follow the surface (11). Moreover, there are (d 2d)/2 linearly independent moving quadrics of degree d -1 that follow the surface and are not obtained from some of the d moving planes by multiplication by x, y, z, w.

The determinant of the (d 2 + d)/2 × (d 2 + d)/2 matrix M d-1 (φ), constructed as in the planar curve case by using the corresponding syzygies, is a power of the implicit polynomial of the surface [START_REF] Cox | On the validity of implicitization by moving quadrics for rational surfaces with no base points[END_REF].

There exist extensions of the method above when the surface has finitely many base points that satisfy certain assumptions, see [START_REF] Busé | Implicitization of surfaces in the projective space in the presence of base points[END_REF].

Turning our attention to the whole syzygy module Syz(f 1 , f 2 , f 3 , f 4 ) of the surface (11), it is not always free, and it is certainly not free when there are no base points. However, if we dehomogenize (11), then the syzygy module is free, of rank 3. Contrary to the case of curves, the elements of the µ-basis are not the syzygies of lowest degree nor they are unique. Moreover we do not have bounds on the degree of the generators. Let µ 3 be the maximum degree. In the case ν ≥ µ 3 , we shall be able to relate ν to the dimension

N ν of Syz(f 1 , f 2 , f 3 , f 4 ) ν .
Let us consider the following map for certain gradings of homogeneous polynomial ring C[s, t, u]:

(C[s, t, u]) 4 ν → C[s, t, u] ν+d : (h 1 , h 2 , h 3 , h 4 ) → 4 i=1 h i f i .
Its kernel is precisely Syz(f 1 , f 2 , f 3 , f 4 ) ν . For ν ≥ ν 0 = 2(d -1), the map is of full rank. Actually, it may be of full rank for lower ν but in constructing implicitization matrices, we are interested in the critical degree. Given a map of full rank, we compute dim Syz(f 1 , f 2 , f 3 , f 4 ) ν as the map's nullity:

N ν = 4 ν + 2 2 - ν + d + 2 2 ,
which is clearly always an integer. This establishes the following.

Lemma 4. For the degree

ν of Syz(f 1 , f 2 , f 3 , f 4 ) ν , ν ≥ 2(d -1) implies dim Syz(f 1 , f 2 , f 3 , f 4 ) ν = 3ν 2 2 -ν(d - 9 2 ) - d(d + 3) 2 + 3.

Orthogonal matrix representations

If we allow for orthogonal matrices and assume that the base points are local complete intersections, then we can restrict ourselves to linear syzygies and construct a matrix M ν (φ) expressing these syzygies, for which the following holds:

Proposition 5. [BLB10, Bus14] Let us define the following critical degrees: ν 0 = d -1 for planar and space curves, and ν 0 = 2(d -1) for triangular surfaces. Then, for all ν ≥ ν 0 , matrix M ν (φ) constructed by the respective linear syzygies satisfies the following property: for any point p ∈ P 2 in the case of planar curves, or p ∈ P 3 otherwise, the rank of M ν (φ) evaluated at p drops if and only if p belongs to the algebraic closure of Im(φ).

We may dehomogenize and obtain the equivalent property, that a point (X, Y ) ∈ C 2 belongs to C if and only if the rank of M ν (X, Y ) drops; the latter denotes the matrix in the non-homogeneous setting.

Syzygies of curves

This section describes how to interpolate the basis of the graded syzygy module of a given degree ν for the case of planar curves, space curves and triangular surfaces. These syzygies can be used to build a matrix as already described.

Planar curves

We describe the method for computing a basis of the linear syzygy module of degree ν of a rational planar curve given by the (unknown) parameterization (3). The dehomogenization of φ gives the rational planar curve C parameterized by

φ : C 1 → C 2 : t → X(t) = f 1 (t) f 3 (t) , Y (t) = f 2 (t) f 3 (t) , ( 14
)
where φ is not known. The input is a set of triplets of the form

(τ 1 ; X 1 , Y 1 ), (τ 2 ; X 2 , Y 2 ), . . . such that φ(τ k ) = (X k , Y k )
, for a range of k ≥ 1 to be defined below. These triplets are assumed sufficiently generic, in particular they may be sampled following the scenarios described in Section 1, e.g. when φ is an arc-length parameterization and the triplets are sampled by a scanner following C.

Our goal is to design an algorithm for computing an implicit matrix representation of the curve C, described by this parametric set of points. The algorithm shall try different degrees ν ≥ 0, and shall compute a C-basis for Syz(X, Y, 1) ν : since the rational functions of X(t), Y (t) are not explicitly known, we compute the basis in the following manner. Consider the moving line h 1 X + h 2 Y + h 3 = 0. The expanded form of each h i is

h i = ν δ=0 h i,δ t δ ∈ C[t], i = 1, 2, 3, (15)
where the h i,δ are (unknown) coefficients. Hence, we can rewrite the moving line as

ν δ=0 t δ Xh 1,δ + ν δ=0 t δ Y h 2,δ + ν δ=0 t δ h 3,δ = 0. ( 16
)
Such equations are going to be used to determine the 3(ν + 1) unknown coefficients h i,δ by interpolation at the sampled triplets. For this, we define a 3(ν + 1) × 3(ν + 1) matrix H whose rows are indexed by evaluations t = τ k , for k = 1, . . . , 3(ν + 1), and each row expresses equation ( 16) as follows:

[X k , τ k X k , . . . , τ ν k X k , Y k , τ k Y k , . . . , τ ν k Y k , 1, τ k , . . . , τ ν k ] .
Clearly, the vector of coefficients [h 1,0 , h 1,1 , . . . , h 3,ν ] corresponding to any element of Syz(X, Y, 1) ν lies in the kernel of matrix H.

We compute a basis of the kernel of matrix H and rewrite the j-th kernel basis vector (h

(j) 1,0 , . . . , h (j) 1,ν , h (j)
2,0 , . . . , h

3,0 , . . . , h

3 ) following equation ( 15). Let h be the kernel dimension of matrix H. We can see that h = N ν under the genericity assumption on the given triplets and the matrix H, because the kernel basis of H corresponds to a C-basis of Syz(f 1 , f 2 , f 3 ) ν . Moreover, since every vector of coefficients of a syzygy lies in the kernel, the vector-space basis of the kernel is a vector-space basis of the syzygy grade, because the latter is a vector space. Then the triplets (h

(j) 1 , h (j) 2 , h (j)
3 ), j = 1, . . . , N ν form a C-basis of Syz(X, Y, 1) ν . In the case h ≥ ν + 1, the C-basis of Syz(X, Y, 1) ν yields the matrix M ν (X, Y ), which offers a matrix representation of the implicit curve C, since ν verifies ν ≥ d -1.

Lemma 2 implies the following. The proof follows easily from the information in Figure 2. Corollary 6. Consider a rational parametric curve C of the form (14). Following the above notation, let d be the homogeneous degree of the (unknown) f i , i = 1, 2, 3, ν ≥ 0 be a fixed degree, specifying a syzygy grading, and h = dim ker(H) be the cardinality of the kernel basis of H. Then,

1. h < ν + 1 if and only if ν < d -1. 2. h = ν + 1 if and only if ν = d -1, then h = d. 3. h > ν + 1 if and only if ν > d -1.
Corollary 6 allows us to compute d by constructing matrix H and comparing h with the selected ν. It is clear that we can also recover the parameterization, but the goal of this work is to obtain robust implicit representations of point cloud models.

Algorithmically, one starts with small ν, say ν = 1. While h < ν + 1, the algorithm doubles ν. If h > ν + 1, we perform binary search to identify the point where h = ν + 1, and h = d.

The first phase, where ν is being doubled, goes up to about 2d, hence needs O(lg d) steps. The binary search takes O(lg d) steps as well, hence the algorithm makes overall O(lg d) corank computations for matrices of dimension up to 2d.

Another possible algorithm uses two values of the syzygy grading, namely ν > ν > 0, and compute the corresponding kernel dimensions h > h ≥ 0. The algorithm terminates when ν ≥ µ 2 -1, then solves h = 2νd + 2 for d. The main step is to compute the slope of the segment defined by the two points, namely

λ = h -h ν -ν ∈ [0, 2].
The algorithm terminates when λ > 1 because this implies ν ≥ µ 2 -1. If λ ≤ 1 the algorithm increases degree ν . This increase happens by setting ν ← ν + 1 then, if λ ≤ 1 again, the algorithm doubles ν . The algorithm requires O(lg µ 2 ) rank computations, which is faster than the previous one.

Rank computation, by means of Gaussian elimination or QR-decomposition, of a m-dimensional matrix has complexity O(m ω ) in the exact setting, where ω < 2.373 is the exponent of matrix multiplication. Clearly, the corank computation for H can be achieved in O(ν ω ) operations, for a given ν. In practice, this is rather of cubic complexity.

The matrices H constructed at various steps are very much related to each other, since the larger ones are obtained by adding columns and rows to a smaller matrix. The new columns and rows correspond, respectively, to higher degree monomials in equation ( 16) and new interpolation points t = τ k . In particular, suppose we have constructed H for some ν, hence of matrix dimension 3(ν + 1), the corresponding matrix H constructed for ν > ν has dimension 3(ν + 1) and the following block structure:

H = H H 12 H 21 H 22 ,
where [H 21 | H 22 ] corresponds to 3(νν) new rows. Suppose the new degree ν = ν + O(1), i.e. the two degrees do not differ significantly, and suppose the corank of H is h = ν -O(1), i.e. it is not significantly smaller than ν. Given a rank revealing decomposition of H, we apply it to the new columns, then compute the rank of H using a total of O(ν 2 ) operations. We thus achieve a speedup of up to one order of magnitude under the current assumptions.

Example 1. Consider the folium of Descartes curve affinely parameterized as:

C = 3t t 3 + 1 , 3t 2 t 3 + 1 ∈ C 2 : t ∈ C (17)
Notice d = 3 for curve C.

Suppose we are given a sample of random points on C for various values of the parameter t, denoted by triplets (τ k ; X k , Y k ), and that we use them to construct the matrix H as described above, with no knowledge of the parametric equation. We try different values of ν:

For ν 1 = 1, the C-basis of Syz(X, Y ) 1 is {(-t, 1, 0)}, that is we are in case 1 of Corollary 6 since N ν 1 < ν 1 + 1. For ν 2 = 2, the computed basis of Syz(X, Y ) 2 is {(-t 2 , t, 0), (-t, 1, 0), (-1/3, -t 2 /3, t)},
that is, case 2 of Corollary 6. This is to be expected since we picked ν 2 = d -1. Any ν ≥ ν 2 is a valid choice to construct an implicit representation matrix M ν (X, Y ).

For ν 2 = 2, the matrix is

M ν 2 (X, Y ) =    -X 0 -Y /3 Y -X 1 0 Y -X/3    , (18)
whose determinant indeed yields implicit equation X 3 + Y 3 -3XY = 0.

Space curves

The method we have described extends naturally to the case of space curves. Assume an unknown parameterization in projective space:

φ : (t 1 : t 2 ) → (f 1 (t 1 : t 2 ), f 2 (t 1 : t 2 ), f 3 (t 1 : t 2 ), f 4 (t 1 : t 2 )), ( 19
)
where d is again defined as the homogeneous degree of the polynomials f i (t 1 , t 2 ), i = 1, . . . , 4. In this case, the critical degree of the syzygies needed for computing the matrix representation of Proposition 5 is d -1, same as for planar curves, meaning ν must be ≥ d -1.

Again, we use moving lines, expressed as follows:

ν δ=0 t δ Xh 1,δ + ν δ=0 t δ Y h 2,δ + ν δ=0 t δ Zh 3,δ + ν δ=0 t δ h 4,δ = 0.
The corresponding equations define matrix H. They contain 4(ν + 1) unknown coefficients h i,δ , hence the dimension of matrix H constructed for some chosen degree ν is dim(H) = 4(ν + 1). Let h be the corank of matrix H.

A corollary of Lemma 3 follows, which shall let us identify the critical degree ν 0 = d -1, see Proposition 5. The proof is straightforward if one considers Figure 1.

Corollary 7. Consider a rational parametric space curve C ⊂ R 3 of the form (19). Let d be the homogeneous degree of the (unknown) f i , i = 1, . . . , 4, ν ≥ 0 be the degree defining the grading of the syzygy, and h = dim ker(H), using the above notation. Then we have:

1. h < 2(ν + 1) if and only if ν < d -1.

2. h = 2(ν + 1) if and only if ν = d -1, then h = 2d.

3. h > 2(ν + 1) if and only if ν > d -1.

We might apply Lemma 3 to establish a similar corollary distinguishing among 3 cases, with middle case h = 2ν + 1. This would have been sufficient for computing d but not enough to build an implicitization matrix from linear syzygies.

Two algorithms are now possible, analogous to those for planar curves in order to identify d and compute the syzygies by interpolation through the parametric point set. Corollary 7 leads to a binary search technique in order to identify the critical degree ν 0 = d -1.

Alternatively, there is an algorithm using two syzygy degrees, namely ν > ν, and computing the slope of the coranks until ν lies in the last segment of the graph in Figure 2. For this algorithm, Lemma 3 implies the following properties for slope

λ = N ν -N ν ν -ν ∈ [1, 3].
First, λ takes an integer value if both ν, ν correspond to the same segment of the polygonal line in Figure 2. Otherwise, we have the following cases:

-λ ∈ (1, 3) iff ν, ν correspond to the first and third segments, -λ ∈ (1, 2) iff ν, ν correspond to the first and second segments, -λ ∈ (2, 3) iff ν, ν correspond to the second and third segments.

Example 2. Consider the Viviani window curve affinely parameterized as:

C = 2t -2t 3 (1 + t 2 ) 2 , 4t 2 (1 + t 2 ) 2 , 1 -t 4 (1 + t 2 ) 2 ∈ C 3 : t ∈ C (20)
The degree of curve C is d = 4.

We are again given a sample of random points on C for various values of the parameter t, denoted by quadruplets (τ k ; X k , Y k , Z k ), which we use them to construct the matrix H.

For ν 1 = 1, the C-basis of Syz(X, Y, Z) 1 is {(-1, -t, t, t), (t, -1, -1, 1)}, that is we are in case 1 of Corollary 7, since N ν 1 < 2(ν 1 + 1). By choosing ν 2 = 3, the computed basis of Syz(X, Y ) 3 consists of 8 elements, that is, case 2 of Corollary 7. This is to be expected since we picked ν 2 = d -1. Thus, any choice of ν such that ν ≥ ν 2 is a valid choice to construct an implicit representation matrix M ν (X, Y, Z).

For ν 2 = 3, the matrix is

M ν 2 (X, Y, Z) =      Z + 1 X/2 0 0 -X/2 0 Y 0 X 1 X/2 0 Z -X/2 2X Y -Y 3X/2 -Y + 1 X -X/2 Z -Y 2X 0 -Y -X/2 -Y -Z + 1 0 -X/2 0 -Y      . (21)

Syzygies of surfaces

This section extends the applicability of our method to surfaces in R 3 without base points.

The theory of syzygies has been fully generalized to certain types of surfaces only, namely tensor product and triangular surfaces [START_REF] Cox | Equations of parametric curves and surfaces via syzygies[END_REF]. In these cases, it is known how many moving planes and moving surfaces, and of which degree, one has to include in order to construct a matrix whose determinant corresponds to the implicit equation. We focus on triangular surfaces because in this case it is easier to obtain the function of the dimension N ν of the graded syzygy module with respect to the degree ν of the syzygies. The method should extend to tensor product surfaces as well, but then N ν is a function of the bi-degree ν = (ν 1 , ν 2 ) of the parameterization.

The input is now a parametric pointset

(τ k , σ k ; X k ), k = 1, 2, . . . , where (τ k , σ k ) ∈ R 2 , and X k ∈ R 3 .
Let us recall triangular surfaces:

φ : P 2 → P 3 : (s : t : u) → (f 1 (s : t : u) : f 2 (s : t : u) : f 3 (s : t : u) : f 4 (s : t : u)) ,
where homogeneous f i (t 1 : t 2 : t 3 ), i = 1, . . . , 4 has degree d. To construct the matrix representation, one has to include d moving planes of degree d -1 and (d 2d)/2 moving quadrics of degree d -1, assuming no base points exist.

To avoid interpolating quadrics and to keep the size of the interpolation matrices low, we shall interpolate only linear syzygies and aim at the critical degree ν 0 = 2(d -1) which, by Proposition 5, allows for constructing an implicit matrix representation by employing only linear syzygies.

As before, it is possible to build an interpolation matrix H, for given degree ν, containing the values of the unknown syzygy monomials at the parametric set of points. The matrix kernel yields the polynomials in the basis of the syzygy grading of degree ν. For a sufficiently generic point sample, the matrix corank h equals the dimension N ν .

Using Lemma 4, namely the quadratic formula

N ν = 3ν 2 2 -ν(d-9 2 )-d(d+3)
2 +3, we can design an algorithm for computing d and interpolate the syzygies beyond the critical degree ν 0 , required for the implicitization matrix of Proposition 5.

The algorithm uses three positive degree values 0 < ν 1 < ν 2 < ν 3 , and computes the 3 respective dimensions N i , i = 1, 2, 3. Then, it checks whether it is possible to fit the 3 pairs (ν i , N i ) on the parabolic formula of N ν as function of ν. If this is possible, we are certain that all 3 values ν i are such that the quadratic formula for N ν holds. Even if N ν as a function of ν is expected to be piecewise with most pieces still known, it is impossible that these 3 points fit another piece, since all pieces are of degree at most 2. Therefore, we can compute d and interpolate the syzygies needed for the implicitization matrix.

Example 3. Consider the canonical Steiner surface affinely parameterized as:

S = 2st s 2 + t 2 + 1 , 2t s 2 + t 2 + 1 , 2s s 2 + t 2 + 1 ∈ C 3 : t, s ∈ C (22)
The degree of the surface S is d = 2.

Given random points on S for various values of the parameters t, s, denoted as 5-tuples of the form (τ k , σ k ; X k , Y k , Z k ), we construct the matrix H.

For ν 1 = 1, the C-basis of Syz(X, Y, Z) 1 is {(-1, 0, t, 0), (-1, s, 0, 0)}, that is we have h = 2. Since we have shown that for ν = d -1 we have h = d linear syzygies, we have successfully computed the degree of the surface, i.e. d = ν 1 + 1 = 2. Thus, any choice of ν such that ν ≥ 2(d -1) = 2 is a valid choice to construct an implicit representation matrix M ν (X, Y, Z).

For ν 2 = 2, the matrix is

M ν 2 (X, Y, Z) =           0 Z 0 0 0 0 0 -X/2 -Y /2 0 0 Z 0 Y 0 0 1 0 -Z/2 0 0 0 0 Y 0 -X/2 0 -X/2 -X 0 Z -X 0 0 0 1 1 0 -X 0 0 -X Y 0 -X/2 -Z/2 0 0 -X 0 0 -X -X/2 -Y /2           . ( 23
)

Implementation and experiments

We experimented using different curves and surfaces of different degrees including the curves we use as examples. All experiments were implemented in Maple 18. The experiments were executed as follows. We start by a given rational parameterization of either a curve or a surface that has no basepoints. That is, we are given a set of 3 or 4 polynomials pols, that is the parameterization of the geometric object in projective space. Then, for random values in the parametric domain we sample the corresponding points on the curve or the surface. We use this parametric pointset for our computations. The parameterization is not used explicitly in our computations apart from verifying the results of our method.

We use this parametric pointset to construct matrix H for a given degree ν, as described in the previous sections. After computing its kernel, we obtain the syzygies that form a basis of the syzygies of degree ν. An example use of our implementation is the command > syzygiesd(pols, t, 3), which returns a basis for the syzygies of degree 3 of the polynomials in pols, whose parameter is t. For different values for the degree ν we look at the number of syzygies we obtain, i.e. the dimension of syzygies of degree ν, and verify their relation to the degree of the parameterization.

The implementation along with the examples included in this paper can be made available upon demand from the authors.

Conclusion and future work

We provide a method for computing a matrix representation of a rational planar or space curve, when we are only given a sufficiently large set of points on the object sampled in such a way that the value of the parameter is known. The algorithm holds for curves without base points in an ambient space of arbitrary dimension, as well as for rational surfaces defined over a triangular patch, again without base points.

One obvious generalization is tensor-product surfaces of bi-degree (d 1 , d 2 ), with parameterization φ : P 1 × P 1 → P 3 : t = (s : t; u : v) → (f 1 (t), . . . , f 4 (t)),

where every f i is bi-homogeneous of degree d 1 in (s : t) ∈ P 1 and degree d 2 in (u : v) ∈ P 1 . In this case, the lack of tight bounds on the degree of the basis of the syzygy module implies that only an approximation to the implicit representation may be obtained.

Future work should involve numerical experiments for interpolating a matrix representation: in this scenario, results are approximate and we wish to quantify the quality of the approximate implicit matrix representation using numerical rank computations. A similar aspect is to consider that noise corrupts the sampling: an estimate of the necessary degree of the syzygies may be obtained in order to interpolate them, thus constructing a matrix approximating the implicit object.

  Then, M ν (φ) is a square d × d implicitization matrix, for ν = d -1. Proposition 1. [CLO05, Sec.6.4],[Cox01, Thm.2.2] When the plane curve has no base points and with the notation above, N d-1 = d and det(M d-1 (φ)) = c • F deg(φ) , where c ∈ C * , F is the implicit polynomial of the curve C, and deg(φ) is the number of pre-images of a generic point on C.

Acknowledgments

The first and third authors are members of "AROMATH", a joint team between INRIA Sophia-Antipolis (France) and National Kapodistrian University of Athens. All authors are partially supported by H2020 Marie Sklodowska-Curie Initial Training Network "ARCADES: Algebraic Representations in Computer-Aided Design for Complex Shapes", 2016-2019. We also thank Laurent Busé for useful discussions during the preparation of this manuscript.