{'fraction_non_alphanumeric': 0.040214098161987476, 'fraction_numerical': 0.011512825691779438, 'mean_word_length': 4.374036687289699, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sampling of Planar Curves: Theory and Fast Algorithms

We introduce a continuous domain framework for the recovery of a planar curve from a few samples. We model the curve as the zero level set of a trigonometric polynomial. We show that the exponential feature maps of the points on the curve lie on a low-dimensional subspace. We show that the null-space vector of the feature matrix can be used to uniquely identify the curve, given a sufficient number of samples. The worst-case theoretical guarantees show that the number of samples required for unique recovery depends on the bandwidth of the underlying trigonometric polynomial, which is a measure of the complexity of the curve. We introduce an iterative algorithm that relies on the low-rank property of the feature maps to recover the curves when the samples are noisy or when the true bandwidth of the curve is unknown. We also demonstrate the preliminary utility of the proposed curve representation in the context of image segmentation.


I. INTRODUCTION
The recovery of a curve from finite number of unorganized and noisy points is an important problem, with applications to computer vision [1], [2] and image processing [3]- [5]. This problem is fundamentally ill-posed because one can find infinite number of curves that pass through the measured points. In addition, the reconstruction is challenging due to curve topology (e.g. the shape of the curve), and the variation of topology with noise. Popular approaches for shape representation with arbitrary topology include (a) explicit representations using a mesh or graphs [1], [6], and (b) implicit level-set representations [7]- [9]. In the first scheme, the shape is constructed from the noisy points as a graph, where the nodes corresponding to adjacent data points are connected. In the second approach, level set functions are constructed from the points [10]. Several methods were introduced to account for noisy data, including spectral graph theory, Laplacian/curvature flow [11], [12]. All of these methods suffer from the inherent parametrization of the curve, which often depends on the sampling density.
The main focus of this paper is to introduce a unified continuous domain theory for the recovery of a planar curve. We assume that the points live on a curve, which is the zero level set of a bandlimited function ψ. This property enables us to express ψ as a finite linear combination of complex exponentials, where the weights are specified by the vector c. The bandwidth of ψ denoted by Λ is a measure of the complexity of the curve [13]. We show that when ψ is irreducible (c.f. Definition 1), the curve consists of a single closed connected component, which we term as an irreducible curve. When ψ can be factorized into multiple irreducible factors, we obtain an union of irreducible curves; each of the irreducible factors correspond to a closed connected curve.
The function ψ vanishes at all points x on the curve; i.e, ψ(x) = 0; ∀x ∈ C. This implies that the weighted linear combination of complex exponential features of the point exp j2πk T x ; k ∈ Λ, weighted by c, will vanish for all points on the curve. In particular, c is the normal vector to the complex exponential features of the points on the curve. We term this property as the annihilation relation, which suggests that the complex exponential maps of the points on the curve lie in a subspace, whose normal vector is c. Thus, the nonlinear exponential mapping transforms the non-linear curve structure to the familiar low-rank or subspace structure, which well-studied in signal processing. When we have a union of irreducible curves, the samples from each one of the irreducible components lie on a subspace; the mapping transforms the complex structure to a union of subspaces structure [14]. The dimension of the subspace spanned by the feature maps is dependent on the bandwidth Λ, and hence the complexity of the curve.
We use the subspace structure of the feature vectors to recover the curve from a few measurements. Specifically, we identify the coefficient vector as the unit norm null space vector of the feature matrix, which is unique up to a scaling with magnitude one. Our worst-case results show that the band-limited function, and hence the curve, can be recovered uniquely, when the number of points exceed a bound that is dependent on the curve complexity. We also show that if the curve is irreducible (single connected component), then one can recover it from any arbitrary sampling pattern. When we have a union of irreducible curves, then the number of samples on each component would depend on the bandwidth of the corresponding irreducible polynomial. We also introduce efficient strategies when the bandwidth of the curve is unknown. Specifically, we show that when the support is overestimated, there exist multiple linearly independent filters that will annihilate the exponential maps; the common zeros, or equivalently the zero level set of the greatest common divisor of the filters, uniquely specifies the curve in this case.
When the support is overestimated, the feature matrix has multiple linearly independent null-space vectors, and hence is low-rank. We note that the Gram matrix of the exponential features correspond to a kernel matrix, which connects the band-limited curve model with widely used non-linear low-rank kernel methods [15]. When the curve samples are noisy, we rely on a nuclear norm minimization formulation to denoise the points. Specifically, we seek to find the denoised curve samples such that their feature vectors form a lowrank matrix. We use an iterative re-weighted algorithm to solve the above optimization problem, which alternates between the estimation of a weight matrix that approximates the null-space and a quadratic sub-problem to recover the data. We note that the iterative algorithm bears strong similarity to Laplacian/curvature flows used in graph denoising, which provides the connection between implicit level-set and explicit graph-based curve representations. One can also derive a graph Laplacian matrix from the weight matrix, which will facilitate the smoothing of signals that live on the nodes of the graph. This graph can be viewed as a discrete mesh approximation to the points that live on the curve. Our experiments show that the Laplacian matrix obtained by solving the proposed optimization algorithm is more representative of the graph structure than classical methods [16], especially when it is estimated from noisy data. This framework reveals links between recent advances in superresolution theory [17]- [19], manifold smoothness based regularization, as well as graph signal processing [20].
This work has connections with Logan's results [21] for the recovery of 1-D band-limited functions from their zero crossings, Fig. 1. Illustration of the annihilation relations in 2D. We assume that the curve is the zero level set of a band-limited function ψ(x), shown by the red function in the top left and the plane slicing the function gives us the level set of the function. The Fourier coefficients of ψ, denoted by c, are support limited in Λ, denoted by the red square on the figure in the bottom right. Each point on the curve satisfies ψ(x i ) = 0. Using the representation of the curve specified by (3), we thus have c T φ Λ (x i ) = 0. Note that φ Λ (x i ) is the exponential feature map of the point x i , whose dimension is specified by the cardinality of the set Λ. This means that the feature map will lift each point in the level set to a Λ dimensional subspace whose normal vector is specified by c, as illustrated by the plane and the red vector c in the top right. Note that if more than one closed curve are presented, each curve will be lifted to a lower dimensional subspace in the feature space, as shown by the two lines in the plane, and the lower dimensional spaces will span the Λ dimensional subspace.
as well as their extensions to 2-D [22]. The main challenge of these works is the extreme sensitivity of the band-limited function to the location of the zero-crossings, when no amplitude information of the signal is used [22]; this has prompted the use of additional information including multi-level contours [22] and multi-scale edges [23]. By contrast, we focus on the recovery of the curve itself, rather than the band-limited function, which is considerably simpler. Specifically, we propose to recover the curve as the zero level set of the sum of squares of all band-limited functions that satisfy the sampling conditions. In addition, unlike [22], our results are also valid for the union of irreducible curves. The proposed work is built upon our prior work on annihilation based super-resolution image recovery [13], [18], [24]- [27] that has similarities to algebraic shape recovery [28] and the recent work by Ongie et al., which considered polynomial kernels [27]. Our main focus is to generalize [27] to shift invariant kernels, which are more widely used in applications. This approach can also be viewed as the generalization of the finite rate of innovation (FRI) theory [29]- [33] in terms of the non-uniform sampling and the knowledge of the samples. We also introduce sampling conditions and algorithms to determine the curve, when the dimension is low. In addition, the iterative algorithm using the kernel trick shows the connections with graph Laplacian based methods used in graph signal processing. The conference version of this work has been published in [34], [35]; this paper provides the proofs of these results, and more elaborate description of details.

A. Parametric level set representation of curves
We model the curve in [0, 1) n , as the zero level set (1) of ψ. We denote the curve specified by the zero level set of ψ by C [ψ]. Note that the level set representation of curves is widely used in segmentation and shape representation [7]. Several authors have proposed to represent ψ as a linear combination of basis functions ϕ k (x) [8], [9]: A popular choice is the shift invariant representation [36] using compactly supported basis functions such as B-splines [8], [9] or radial basis functions [10]. In this case, we have ϕ k (x) = ϕ(x − k), where ϕ(x) is a B-spline or Gaussian function. Here, |Λ| denotes the number of basis functions, which is equivalent to the number of grid points in the context of shift-invariant representations.

B. Band-limited level set representation of curves
We now consider a special case of the parametric level set curve representation, in which the basis functions are chosen as ϕ k = exp(j 2πk T x); k ∈ Λ, j = √ −1. Then Here, the Fourier coefficients c k are supported on a rectangular grid Λ ⊂ Z n . We use Hermitian symmetry property to create the curves in our experiements. But in the proofs and estimation, the Hermitian symmetry property is not utilized. Here, BW(ϕ) = |Λ| is the bandwidth of the curve C[ψ]; |Λ| denotes the cardinality of the set. See Fig. 1 for an illustration. Note that for each index k, the basis functions satisfy |ϕ k (x)| = 1; ∀x. The non-compact nature of the basis functions unlike [8], [9] makes the analysis of the problem easy. We note that the band-limited level set model can represent closed curves with arbitrary complexity. See Fig. 2 in [18] for examples. The complexity of the curve is dependent on the cardinality of the set Λ, which we denote by |Λ|, which is the number of free parameters in the curve representation. We now illustrate the geometric properties of bandlimited curves.

C. Properties of bandlimited curves
We now briefly introduce the relation between band-limited functions and complex polynomials [18] in 2-D. We use the change of variables: (x1, x2) → (z1, z2) = (exp(j2πx1), exp(j2πx2)), which gives a one-to-one mapping [18] P from [0, 1) 2 to the unit torus T 2 . This mapping allows us to rewrite (3) as the complex polynomial: Since P is one-to-one, the space of bandlimited functions ψ(x) in (3) and the space of complex polynomials living on T 2 are isomorphic. We can analyze the properties of P[ψ] to deduce the properties of the trigonometric polynomial (3).
1) Irreducible curves: We now consider the simplest possible curve, which consists of a single connected component.
Definition 1 (Irreducible trigonometric polynomials). A polynomial η(z) with complex coefficients is irreducible, if it cannot be expressed as the product of two or more non-constant polynomials with complex coefficients [37]. The trigonometric polynomial η(x) is termed as irreducible, if the polynomial specified by P[η] is irreducible.
Definition 2 (Irreducible curve). A trigonometric curve C is irreducible if it is the zero level set of an irreducible trigonometric polynomial.
In most cases, the zero level sets of irreducible trigonometric polynomials have only one component, except for some pathological cases. We consider the irreducible curves with only one connected component in the rest of the paper. Note that by definition, the curves are periodic in [0, 1) 2 . We observe from the segmentation examples in Fig. 9 that this restriction does not pose a problem in practical applications; some additional boundary curves are added, which can be ignored in practice. If this poses a problem, the image may be padded by additional zero regions. While the representation cannot represent open curves in the strict sense, an open curve can be approximated with arbitrary accuracy by a bounded closed curve [38].
2) Union of irreducible curves: The above curves are often too simple to represent real-world structures. One can represent curves with multiple closed & connected components as the zero level set of a bandlimited function with multiple irreducible factors ψ = η1 · η2 . . . ηn. In this case, the resulting curve is the union of the corresponding irreducible curves Note that the coefficient vector of the union of irreducible curves ψ is obtained as the convolution between the coefficient vectors of the irreducible components ηi: where di F ↔ ηi. This implies that the bandwidth of ψ is greater than that of the components.
3) Minimal Polynomial: Note that if ψ(x) defines a curve, then η(x) = exp(j2πk T x) · ψ(x) will also define the same curve for k ∈ Z 2 since exp(j2πk T x) never vanishes in [0, 1) 2 . But the bandwidth of η is higher than that of ψ. We term the polynomial with the smallest bandwidth to be the minimal polynomial of a curve, which is the object we are interested in in this paper. The following result about minimal polynomial of a curve guarantees the existence of the minimal polynomial. We note that in the 1-D setting (n = 1), there is a one-to-one correspondence between the number of points on the zero set of ψ and the bandwidth |Λ| [18], [28]; this relation enabled the use of |Λ| in [13], [18] as a surrogate for sparsity in signal recovery. In higher dimensions, |Λ| can still serve as a complexity measure. Note that the zero-set can consist of isolated Diracs in 2D, which implies that BW(ψ) can still serve as the surrogate for sparsity. However, we emphasize that (3) can provide a significantly richer representation, even when the signal is not isolated and consists of points on a curve.

A. Annihilation relations for points on the curve
Consider an arbitrary point x on the curve specified by (1) and (2). By definition, we have ψ(x) = 0, which translates to: . . .
Note that φΛ : R n → C |Λ| is a non-linear mapping or lifting of a point x to a high dimensional space, whose dimension is given by the cardinality of the set Λ, denoted by |Λ|. Note that this non-linear lifting strategy is similar to feature maps used in kernel methods. We hence term φΛ(x) as the feature map of the point x. Note that every point on the curve satisfies (7), which we term as the annihilation relation.
Let us now consider a set of N points on the curve, denoted by x1, · · · , xN . Note that the feature maps of each one of the points satisfy the above annihilation relations, which can be compactly represented as: Here, ΦΛ(X) is the feature matrix of the points and X = [x1 x2 . . . xN ].
Assume that we have a union of irreducible curves as in (5), where the bandwidth of each of the irreducible components C[ηi] is Λi and bandwidth of C[ψ] is Λ. In this case, the |Λi| dimensional lifting ΦΛ i (x) of the samples on C[ηi] will lie on a |Λi| − 1 dimensional subspace. Similarly, the |Λ| dimensional lifting ΦΛ(x) of the samples on the union of irreducible curve C[ψ] will lie on a |Λ|−1 dimensional subspace.

B. Curve recovery from samples
When ΦΛ is rank-deficient by one, the coefficient vector c can be identified as the unique non-zero null-space basis vector of ΦΛ(X). This implies that the features lie in an |Λ| − 1 dimensional subspace, whose normal is specified by c. This annihilation relation is illustrated in Fig 1, in the context of band-limited curves considered in the next subsection. We will show that there exists a unique null-space basis vector when complex exponential basis functions are chosen as in Section II-B.
In practice, the points are often corrupted by noise. In the presence of noise, the null-space conditions are often not satisfied exactly. In this case, we can pose the least square estimation of the coefficients from the noisy data points {xi} N i=1 as the minimization of the criterion: where QΛ = N i=1 φΛ(xi)φΛ(xi) T . To eliminate the trivial solution c = 0, we pose the recovery as the constrained optimization scheme: The solution is the eigenvector corresponding to the minimum eigenvalue of QΛ. Note that QΛ is nothing but ΦΛ(X)Φ T Λ (X). Thus we just need to use the singular value decomposition of Φ T Λ (X) to obtain the desired solution.

C. Irreducible band-limited planar curve: sampling theorem
We now focus on the problem of the recovery of the curve (1), given a few points {xi ∈ R 2 ; i = 1, · · · , N } on the curve 1 . Let us take the band-limited curve representation as (3) for the rest of the section to derive our sampling conditions. We now determine the sampling conditions for the perfect recovery of the curve ψ(x) = ↔ ψ is support limited to a 3 × 3 region, shown in (a). The theorem guarantees the perfect recovery will happen if we have no less than (k 1 + k 2 ) 2 = 36 samples. We first randomly chose 36 samples on the curve. Then from these 36 randomly chosen samples, we obtained (b), which gives us perfect recovery of the original curve. Furthermore, we mentioned that we do not require any constraint on the distribution of samples on the curve. In (c), we randomly chose 36 samples from the left half part of the curve and we got perfect recovery as well. In (d), 36 samples are randomly chosen from the right half of the curve. From (d), we saw that perfect recovery of the whole curve was also obtained. For each case, the average time required for the recovery is about 1.2 second. 0 using (10), when the curve is specified by (3). In this case, the annihilation relation (7) is satisfied with the feature maps defined as We also assume that Λ is a rectangular neighborhood in Z 2 of size k1 × k2. We first review some results from algebraic geometry.
There is a one-to-one correspondence between trigonometric polynomials and complex polynomials, as shown in Section II-C. We use the extension of Bézout's inequality for trigonometric polynomials, which bounds the number of solutions of the system µ(x) = η(x) = 0 that do not have any common factors.
Lemma 5 (Bézout's inequality for band-limited polynomials). Let µ(x) and η(x) be two band-limited polynomials, whose Fourier coefficients are support limited to k1 × k2 and l1 × l2, respectively. If µ and η have no common factor, then the system of equations has a maximum of (k1 The proof of Lemma 5 is given in Appendix VII-A. We use this property to derive our main results. We first focus on the case where ψ is an irreducible band-limited function. be N distinct points on the zero level set of an irreducible band-limited function ψ(x), x ∈ R 2 , whose Fourier coefficients are restricted to a rectangular region Λ with size k1 ×k2. Then the curve ψ(x) = 0 can be uniquely recovered by (8), when: The proof is provided in Appendix VII-B. Note that the sampling condition for a single irreducible curve does not specify any constraint on the distribution of points on the curve; any set of N > (k1 +k2) 2 points are sufficient for the recovery of the curve. This proposition is illustrated in Fig. 2, which shows that the recovery is guaranteed irrespective of the distribution of samples. This property is similar to well-known results in non-uniform sampling of band-limited signals [39], where the recovery is guaranteed under weak conditions on the nonuniform grid and the average sampling rate exceeding Nyquist rate.
We compare this setting with the sampling conditions for the recovery of a piecewise constant image, whose gradients vanish on the zero level set of a band-limited function [18]. The minimum number of Fourier measurements required to recover the function there is |3Λ|. When k1 = k2 = K, then 9K 2 complex Fourier samples are required, which is far more than 4K 2 real samples required for the recovery of the curve in our setting. Note that the constant values within the regions bounded by the curves also need to be recovered in [18], which explains the higher sampling requirement. We note that the above bounds are looser than the ones in [32], which are based on the number of available equations; they assume the chances of the equations being linear dependent is unlikely [32]. Unlike the high-probability results in [32] that our bounds are worst-case guarantees, which will hold irrespective of the sampling geometry. We note from the experiments in Fig. 5 that recovery succeeds in most cases whenever the number of samples exceed |Λ| − 1 = k1k2 − 1, which is the number of degrees of freedom in representing the curve.

D. Union of irreducible curves: sampling theorem
We now generalize the previous result to the setting where the composite curve is a union of multiple irreducible curves. Equivalently, the level set function is the product of multiple irreducible band-limited functions. We have the following result for this general case: , where the bandwidth of the j th factor is given by k1,j × k2,j. The curve ψ(x) = 0 can be uniquely recovered by (8), when each of the irreducible curves are sampled with The total number of samples needed for unique recovery is specified by which is bounded above by (k1 + k2)(k1 + k2 + 2(J − 1)).
We note that the upper bound can be approximated as (k1 + k2) 2 for small values of J, which is the upper bound in Proposition 6. The above result is proved in Appendix VII-C. Note that unlike the case considered in Section III-A, an arbitrary set of N samples cannot guarantee the perfect recovery. Each of the J irreducible curves C[ηj] need to be sampled proportional to their complexity, specified by deg(ηj) to guarantee perfect recovery.
We demonstrate the above proposition in Fig. 3. We consider a curve C[ψ], where c ψ F ↔ ψ is support limited to a 5 × 5 region. We note that there are three connected components in the above curve. We consider the recovery from different number of samples of C[ψ] in the middle row, sampled randomly. The random strategy ensures that the samples are distributed to the factors, roughly satisfying the conditions in Proposition 7. Note that the theory guarantees recovery, when the number of samples exceeds around (k1 + k2) 2 = 100 samples. We observe good recovery of the curve around 50 samples; note that our results are worst-case guarantees, and in practice fewer samples are sufficient for good recovery of most curves. We further study the distribution of the points in Fig. 4. The experiments demonstrate that each of the curves need to be sampled with a number  Fig 3 (b), which is given by the union of two irreducible curves with bandwidth 3 × 3. So the bandwidth of C[ψ] is 5 × 5. According to Proposition 7, we will need to have around 100 samples to recover C[ψ] and each of the two irreducible curves need to satisfy with the sampling condition. As we noted in Fig 3, our results are worst-case guarantees. We observe that when we have 50 points and those points are uniformly sampled on the two irreducible curves, we can successfully recover the whole curve, as shown in (a). Now, if we put most of the samples on one of the irreducible curves, we cannot fully recover the curve, as illustrated in (b) and (c). This implies that the sampling condition on each of the irreducible factors is necessary in Proposition 7. proportional to the bandwidth of the curves as in (b). When the points are non-uniformly distributed as in (c) or (d), the recovery fails. We further studied the above proposition in Fig. 5. We considered several random curves, each with different bandwidth and considered their recovery from different number of samples. The sampling locations were picked at random. The colors indicate the average reconstruction error between the actual curve and the reconstructed curves. This reconstruction error is computed as the sum of distances between each point on one curve and the closest point to it on the other curve. We have also plotted the upper bound (k1 + k2) 2 in red, while the number of unknowns in the curve representation k1 k2 is plotted in blue. We note that the curve can be recovered accurately when the number of samples exceed the upper bound. We also note that in general, good recovery can be obtained for most curves, when the number of samples exceed k1 k2.
E. Curve recovery with unknown Fourier support Propositions 6 and 7 assume that the true support of the Fourier coefficients of ψ, specified by Λ is known, in addition to the points {xi} N i=1 . However, typically only the points will be known and the filter support will be unknown. We now consider the case where the filter support is over-estimated as Γ ⊃ Λ. We focus on the recovery of the coefficients from the annihilation relation The following result shows that the above matrix will have multiple linearly independent null-space vectors. However, if the curves are sampled as described below, the corresponding band-limited functions satisfy some desirable properties that facilitate the recovery of the curves.
all of the above functions, or equivalently the right nullspace vectors cµ F ↔ µ of ΦΓ, will be of the form: where η(x) is an arbitrary function such that supp(cµ) = Γ.
Note that the minimal function ψ(x) is a special case of (18), with η = 1. The above result is proved in Appendix VII-D. Since ψ(x) is the common factor of all the annihilating functions, all of them will satisfy µ(x) = 0, for any point on the original curve as well as the sampling locations. This also implies that ψ(x) is a Fig. 5. Effect of number of sampled points on perfect reconstruction. We randomly generated several curves with different bandwidth and number of sampled points, and recovered the curves from these samples. The success of reconstruction of the curves averaged over several trials are shown in the above phase transition plot, as a function of bandwidth and number of sampled entries. The color indicates the frequency of success; the color black indicates that the true curve cannot be recovered in any of the experiments, while the color white represents that the true curve is recovered in all the experiments. It is seen that perfect recovery occurs whenever we have ≥ (k 1 + k 2 ) 2 samples, as indicated by our worst-case guarantees. However, we note that good recovery is observed whenever the number of samples exceed the degrees of freedom k 1 · k 2 . We consider the recovery of the curve C[ψ] as specified by Fig 3 (b), assuming unknown bandwidth. We overestimate the support Γ as 11x11, while the original support of c ψ F ↔ ψ is 5 × 5. According to Propositions 8& 9, when the number of samples exceed (k 1 + k 2 )(l 1 + l 2 ) = 220, the matrix is low-rank. The first row shows the results by using 220 samples. We display the Fourier transforms of the three null-space functions of Φ Γ in (a), (b) and (c). This approach of visualizing the null-space functions is similar to the approaches in [32], [40]. All of these functions are zero on the C[ψ], in addition to possessing several other zeros. The sum of squares function, denoted by (22) is shown on the right column, captures the common zeros, which specifies the curve C[ψ]. We use the SOS function as a surrogate for the greatest common divisor of the null-space functions. Note that the bound in Proposition 8 is also a worst-case guarantee. In the second row, the curve C[ψ] was sampled on 100 random sampling locations, denoted by the red crosses. We see that the curve can be recovered well using just 100 samples. The computational time used to specify the curve using SOS function in this experiment is about 1.6 second . common divisor of the above functions µ(x). In fact, ψ(x) is the greatest common divisor as we will show it in the next paragraph. We now characterize the number of linearly independent annihilation functions, or equivalently the size of the right null space of ΦΓ. with equality if the sampling conditions of Proposition 8 are satisfied. Here, This set is illustrated in Fig 5(a) from [18] along with Γ and Λ. The inequality of this result is same as the inequality of Proposition 5.1 in [18]. Based on the inequality, we can then obtain the second part (the equality) of the result, which provides us a means to compute the original curve, even when the original bandwidth/support Λ is unknown. Specifically, Proposition 9 shows that ΦΓ(X) has |Γ : Λ| null-space vectors, each of which satisfies (18). Besides, from the proof of Proposition 9, we can see that any polynomial of the form is a null-space vector of ΦΓ(X). Note that the exponentials exp(j2πl T x), ∀l ∈ Γ : Λ are linearly independent, and hence the set {θ l ; l ∈ Γ : Λ} spans the null space of ΦΓ(X). Since exp(j2πl T x) does not vanish in the domain, the only common zeros of {θ l ; l ∈ Γ : Λ} will be the zeros of the minimal polynomial ψ(x), meaning that ψ(x) is the greatest common divisor of the functions that span the null-space of ΦΓ(X). Therefore, the common zeros of these functions, or equivalently the zeros of the greatest common divisor, will specify the curve. A cheaper alternative to evaluating the greatest common divisor is to evaluate the sum of squares polynomial, specified by: which will vanish only on points satisfying ψ(x) = 0. Here Q = |Γ| − r is the dimension of the right null-space of ΦΓ. Since c ψ F ↔ ψ is a valid right null-space vector of ΦΓ that only vanishes on the true curve, the sum of squares function γ specified in (22) will only vanish on the true curve. Thus, if the total number of points sampled are N = J j=1 Nj > (l1 + l2)(k1 + k2 + 2(J − 1)), and are arranged as (17), then the curve can be uniquely recovered.
We demonstrate the above result in Fig. 6. We considered the sampling of the same curve illustrated in Fig. 3, with the exception that we over-estimated the support to be 11 × 11 as opposed to the true support of 5 × 5. We considered 220 random samples, which satisfies the sampling conditions in Proposition 8. We show three of the annihilating functions in the first three columns of Fig. 6. We note that all of these functions are valid annihilating functions, but possess additional zeros. By contrast, the sum of square polynomial shown on the right uniquely specifies the curve.

IV. APPLICATIONS
The above section provided details on curve representation as well as provided guarantees on when the curve can be recovered from finite number of samples. We will now demonstrate the utility of the theory in some preliminary representative applications. We believe that the representation has the potential to improve upon the state of the art in several areas. However, these developments and rigorous validations are beyond the scope of this preliminary work. The experiments in this section are run on a laptop with Intel Core i7-9750H CPU.

A. Application in segmentation
The Mumford Shah functional is a popular formulation for segmenting objects into piecewise constant regions. It approximates an image f by a piecewise constant function in the 2 sense, where Ω k , k = 1, .., K are the regions and a k are the constants and χ represents the characteristic function on the set. The bounded curve is denoted by ∂Ω = ∂Ω1 ∪∂Ω2 ∪· · ·∪∂ΩK . Different penalties, including the length of ∂Ω or its smoothness are imposed to regularize the optimization problem. We propose to represent ∂Ω as the zero level set of a band-limited function ψ specified by (3) as in [18]. In this case, the piecewise constant function satisfies ∇f * c = 0, which can be expressed in the matrix form as where T is a block Toeplitz 2-D convolution matrix andc is the matrix version of vector c. When the bandwidth is over-estimated, T ∇f has multiple linearly independent null-space vectors and hence the matrix is low-rank. Note that the rank of the matrix can be considered as a surrogate for the complexity of the curve ∂Ω. We hence formulate the segmentation task as the low-rank optimization problem, analogous to [40].
where h is the original image. Note that as λ → ∞, T ∇f approaches a rank r matrix. Once f * is obtained, the sum of square function of the null space of T ∇f will specify the curve and f * is the piecewise constant approximation. We use an alternating minimization strategy as reported in [40] to solve the above optimization scheme.

B. Recovery of noisy point clouds
We now consider the case where we have noisy measurements of points lying on the curve ψ(x) = 0, where ψ is represented as a linear combination of basis functions as in (7). When the measurements are noisy, we propose to denoise them using the low-rank property of the features discussed in Proposition 9.
1) Relation to kernel methods: Proposition 9 indicates that we can solve inverse problems by enforcing a low-rank constraint on the feature matrix. However, the size of the feature matrix grows with the dimensionality of the ambient space n as well as the size of the overestimated filter support Γ. Thus, in practice, it might be infeasible to form the feature matrix. However, the Gram matrix given by is of size N × N , where N is the number of points. Note that the complexity of an algorithm that depends on the Gram matrix is independent of the dimension of the ambient space and the chosen filter support Γ. We note that this approach is similar to the "kerneltrick" used in various machine learning applications [15]. Under our assumed model, the entries of this Gram matrix are given by When Γ is a centered cube in R n , κΓ(r) is a Dirichlet function. Note that the entries of KΓ can be evaluated without explicitly evaluating the feature matrix. The width of the Dirichlet function is dependent on the Fourier support Γ. The kernel matrix satisfies rank(KΓ) ≤ r, where r is given by (19). The relationship (26) implies that the kernel matrix is low-rank, which is a property that is widely used in kernel low-rank methods. A popular practical choice for kernel is Gaussian functions, or equivalently periodized Gaussian functions when the domain is restricted to [0, 1) 2 . Note that the Gaussian kernel is qualitatively similar to the Dirichlet kernel considered above, where the width of the Gaussian is a parameter similar to the size of the support Λ. We note that the Gaussian kernel function is less oscillatory and is isotropic, which makes it more attractive than Dirichlet 2 kernels in applications. The Gaussian kernel correspond to feature maps of the form as Γ → Z n . Note that this is the Fourier transform of a shifted Gaussian; this setting corresponds to the case where the level set function is being expressed as a shift invariant linear combination of Gaussian functions on a very fine grid. Since the Gaussian kernel matrix KΓ is theoretically full rank, the theoretical analysis in the previous section is not directly applicablele to Gaussian kernels in the strict sense. However, we observe that the Fourier series coefficients of a Gaussian function may be approximated to be zero outside |k| < 3/πσ, which translates to |Λ| ≈ 6 πσ n . This implies that the Gaussian kernel matrix can be safely approximated to be low-rank, when σ is sufficiently high; this corresponds to a more localized kernel in space.
In the next subsection, we will describe how to use proposition 9, without explicitly forming the feature matrix ΦΓ(X), and computing only it's Gram matrix KΓ instead.
2) Denoising of point clouds using nuclear norm minimization: With the addition of noise, the points deviate from the zero level set of ψ. A high bandwidth potential function is needed to represent the noisy curve. Let the noisy measurements of the matrix X be given by Y. We propose to use the nuclear norm of the feature matrix as a regularizer in the recovery of the points from noisy measurements: We note that this approach has conceptual similarities to Cadzow denoising that is widely used in FRI methods [31], [32]. We use an IRLS approach [41], [42] where the nuclear norm is approximated as: where P = [K(X) + γI] − 1 2 and K is the Gaussian kernel. Here, γ is a small constant added to ensure that the inverse is well-defined.
Note from (31) that updating X involves the solution of a nonlinear system of equations. We propose to linearize the gradient of the cost function in (31) with respect to X.
Linearizing the gradient with respect to Xi, we obtain: where w In matrix form, the gradient can be rewritten as Here, D (m) is a diagonal matrix with elements defined as D ij . This results in the following equivalent optimization problem for the estimation of X at the n th iteration, which can be solved analytically: Thus, we alternate between the estimation of X (m) and L (m) till convergence to solve (29). We note that this optimization algorithm is non-convex and does not come with any convergence guarantees. However, with reasonable initialization (e.g. X = Y), the algorithm yielded good results in practice.

V. COMPARISONS WITH STATE OF THE ART
Now, we compare our planar curves recovery results with an algorithm that relies on level set evolution. Specifically, we re-engineer the level-set based method for curve recovery termed as "distance regularized level set evolution" (DRLSE), which was introduced in [43] for image segmentation. DRLSE poses the image segmentation as the minimization of the cost function where φ is the level set function. Here, Rp(φ) is a level-set regularization term which maintains the level-set function φ as a signed distance function. We choose the function to be (16) of [43]. The first and second terms are the weighted length and area of the curve, respectively: which are determined by the choice of edge indicator function g. Length and area minimizing flows are well-studied in the level-set literature, and correspond to curve velocities that are proportional to curvature and constant velocity along the curve normals [44]. The parameter µ for level set regularization term is determined by the time step. Once the time step is chosen, µ is almost determined because of the Courant-Friedrichs-Lewy (CFL) condition. We will re-engineer DRLSE to the curve recovery from samples by choosing the edge indicator function as the distance of the level-set function to points. Since DRLSE was originally designed for image segmentation, we use the edge-based edge indicator function discussed in (40). We call the re-engineered DRLSE algorithm the level-set based algorithm.

A. Curve recovery from samples
In this section, we compare the proposed curve recovery scheme in Section III with the level-set based algorithm. We choose the edge indicator function as the distance of the curve from the samples xi; i = 1, .., N : where d(x, xi) = mini=1,··· ,N {c, x − xi 2 } for all x in the image domain and c is a large constant. We compare the two methods in the context of recovering the edge curve for the Chinese character "Tian" (meaning sky in English) in Fig. 7.
In our method, we chose the bandwidth of the curve as 51×51. For the level-set based algorithm, we choose the parameters as λ = 5, α = 10 and the initialization curve as a square which includes the whole curve. Note that our method do not need any initialization. The two rows show the recovery results by the two different methods from 600 and 1000 samples respectively. The first column shows the samples we choose. The results obtained by using the level-set based algorithm are given in the second column. The numbers of iterations for obtaining (b) and (e) are 1510 and 2260. The third column shows the recovery results by using our method. By comparing (b) and (c), one can see that our method recover the curve successfully from 600 randomly chosen samples. For the level-set based algorithm, the curve is not successfully recovered from those 600 samples. Once we have 1000 samples, we can find that both the two methods succeed in recovering the curve, as shown in (e) and (f). However, the computational time required for our proposed algorithm is less than that of the level-set based algorithm. This example also demonstrates that both the two level-set based methods work well even though the curves have some sharp corners.

B. Segmentation using structured low-rank methods
We demonstrate the preliminary utility of the scheme in Section IV-A in Fig. 9 on two images: the cells image and the church image. The proposed algorithm in (25) is initialized with f = h and iterated until convergence. The parameter λ is set to a high value (e.g. 5×10 9 in our experiments) to enforce the rank constraint. We note that the  9. Illustration of sensitivity of our proposed image segmentation algorithm to the rank. From the segmentation results, we see that when the rank is small, simpler segmentation curves will be obtained. When the rank is chosen to be too high, we will obtain over-segmentation result. Thus, the rank is a good surrogate for the complexity of the curve.
optimization scheme is capable of identifying the cells and the outline of the church, even though no curve initialization was provided.
For the cells segmentation (c) and the church segmentation (f), we choose the rank to be 500 and 1200 respectively. The corresponding computational time for the two segmentations is 60 seconds and 226 seconds. We now compare our segmentation method with the levelset based segmentation method DLRSE, where the edge indicator function is chosen as where Gσ is a Gaussian kernel with a standard deviation σ. We considered the initialization of DLRSE with two possible curves, indicated by the green squares. The parameters in DLRSE were chosen manually to yield the best results, which corresponded to These results show that the proposed scheme is comparable to DLRSE in segmentation performance and can capture sharp features. However, the main benefit is its insensitivity to initialization, compared to DLRSE seen from (d) and (e).

C. Point cloud denoising: comparison with state-of-the-art
We illustrate the utility of the kernel low-rank formulation (KLR) proposed in IV-B2 to denoise 2D points in Fig. 10. In Fig. 10, we also compared our method with another point-set denoising method called Graph Laplacian regularized point cloud denoising (GLR), which was introduced recently in [45]. We choose 3 examples to perform the point sets denoising algorithms and in each example, we add Gaussian noise to the point sets. In the first example, we randomly choose 409 points on the edge set of a rabbit as shown in (a). We use GLR and KLR to denoise the noisy points respectively. For GLR, we set the parameter µ to be 1000 and after 34 iterations, we obtained the denoising result (c). In KLR, we get the denoising result (d) after 80 iterations using about 4.4 seconds. In the second examples, we choose 385 points on the edge set of a plane. We again set the parameter µ to be 1000 and after 31 iterations, we have the denoising result (g). For KLR, we get the result (h) by iterating 80 times using about 4.0 seconds. In the third example, we choose 451 points on the shape of a fish. For GLR, after 34 iterations by setting the parameter µ in the algorithm to be 1000, we obtain the denosing result (k). For our proposed denoising algorithm, we raised the number of iterations to 450 and it takes about 32.7 seconds to obtain the denoising result (l). By comparing the denosing results (c) and (d), (g) and (h), (k) and (l), we can see that both the two methods work for denosing the noisy points. While for GLR, we can see that the some points will get closer along the right curve. To compare the experimental performance mathematically, we introduce an evaluation metric, signal-to-noise ratio (SNR), for point cloud denoising. Suppose the ground-truth and predict point clouds are {xi} N 1 i=1 and {yi} N 2 i=1 . We define the SNR, which is measured in dB by where MSE is the mean-square-error defined as

VI. DISCUSSION AND CONCLUSION
We introduced a continuous domain framework for the recovery of points on a band-limited curve. The proposed bandlimited representation have several desirable geometric properties, which make it an attractive tool in a variety of shape estimation problems. We have introduced novel algorithms with sampling guarantees for the recovery of both irreducible and union of irreducible bandlimited curves from few of their samples; our experiments show that the bounds are sharp. We can see from Fig. 5 that with less samples, perfect recovery will also happen. We will discuss this phenomenon in our feature work. Furthermore, we do not use the property that the coefficients of the polynomial are Hermitian symmetric [32] in the proofs and estimation. We note that the use of the Hermitian symmetry property may result in improved sampling bounds. However, this extension is beyond the scope of the present work.
We also demonstrated the utility of the representation in practical applications including image segmentation and denoising of a point cloud, which can be modeled by a curve. The main benefit of the curve recovery from points as well as image segmentation over the state of the art is the convex formulation, which makes the algorithm insensitive to local minima errors as well as initialization. The segmentation and point cloud denoising experiments show that the proposed scheme can exploit the global structure of the points better than competing methods that rely on local curve properties such as smoothness and curvature, which makes the algorithms less sensitive to non-uniformity of sampling. In our future work, we will extend our results to high dimensional spaces. We note that the proofs based on Bézout's inequality cannot be readily extended to the higher dimensional setting and thus a new method need to be introduced for the extension.

ACKNOWLEDGEMENT
The authors gratefully acknowledge helpful comments from Dr. Greg Ongie.

A. Proof of Lemma 5
We first state the well-known result for complex polynomials, which we extend to the band-limited setting.

B. Proof of Proposition 6
Proof. The Fourier coefficients of ψ(x) is support limited within Λ, which is the minimal support. Let η(x) be another band-limited polynomial, whose Fourier coefficients are support limited within Λ and satisfies η(xi) = 0, for i = 1, . . . , N . When the number of samples satisfy (13), this is only possible if η is a factor of ψ, according to Bézout's inequality. Thus, ψ(x) must be a factor of η(x). Since ψ is irreducible, this implies that it is the unique bandlimited irreducible polynomial satisfying ψ(xi) = 0.

C. Proof of Proposition 7
Proof. The polynomial ψ(x) is represented in terms of its irreducible factors as: ψ(x) = ψ1(x)ψ2(x) . . . ψJ (x) where the bandwidth of ψj(x) is k1,j × k2,j. Let η(x) be another polynomial with bandwidth k1 × k2 satisfying η(xi) = 0, for i = 1, . . . , N . Consider one of the irreducible subcurves {ψj(x) = 0}, that is sampled on Nj points satisfying (14). According to Lemma 5, both ψj and η can be simultaneously zero at these sampling locations only if ψj and η have a common factor. Since ψj is irreducible, this implies that ψj is a factor of η. Repeating this line of reasoning for all factors {ψj}, we conclude that ψ(x) divides η(x). Since both ψ(x) and η(x) have the same bandwidth, the only possibility is that η(x) is a scalar multiple of ψ(x). This implies that the curve ψ(x) = 0 can be uniquely recovered in (14) is satisfied.

D. Proof of Proposition 8
Proof. Following the steps of the proof for Proposition 7, we can conclude that ψ(x) is a factor of µ(x). Since Λ ⊂ Γ, it follows that µ(x) = ψ(x) η(x), where η(x) is some arbitrary function such that µ(x) is band-limited to Γ.

E. Proof of Proposition 9
Proof. Let c be the minimal filter of bandwidth |Λ|, associated with the polynomial ψ(x). We define the following filters supported in Γ for all l ∈ Γ : Λ.
If the sampling conditions of Proposition 8 are satisfied, then all the polynomials corresponding to null-space vectors of ΦΓ are of the form: µ(x) = ψ(x) η(x). Alternatively, in the Fourier domain, the filters are of the form: where d l are the Fourier coefficients of the arbitrary polynomial η(x). Thus, all the null-space filters can be represented in terms of the basis set {c l }. This leads to the relation: rank (ΦΓ(X)) = |Γ| − |Γ : Λ|.