Sampling and Reconstruction of Shapes with Algebraic Boundaries

Sampling and Reconstruction of Shapes with Algebraic Boundaries

Mitra Fatemi, Arash Amini, and Martin Vetterli M. Fatemi and M. Vetterli are with the Department of Computer and Communication Sciences, École Polytechnique Fédéral de Lausanne (e-mail: {mitra.fatemi, martin.vetterli} Amini is with the Department of Electrical Engineering, Sharif University of Technology (e-mail:

We present a sampling theory for a class of binary images with finite rate of innovation (FRI). Every image in our model is the restriction of to the image plane, where denotes the indicator function and is some real bivariate polynomial. This particularly means that the boundaries in the image form a subset of an algebraic curve with the implicit polynomial . We show that the image parameters –i.e., the polynomial coefficients– satisfy a set of linear annihilation equations with the coefficients being the image moments. The inherent sensitivity of the moments to noise makes the reconstruction process numerically unstable and narrows the choice of the sampling kernels to polynomial reproducing kernels. As a remedy to these problems, we replace conventional moments with more stable generalized moments that are adjusted to the given sampling kernel. The benefits are threefold: (1) it relaxes the requirements on the sampling kernels, (2) produces annihilation equations that are robust at numerical precision, and (3) extends the results to images with unbounded boundaries. We further reduce the sensitivity of the reconstruction process to noise by taking into account the sign of the polynomial at certain points, and sequentially enforcing measurement consistency. We consider various numerical experiments to demonstrate the performance of our algorithm in reconstructing binary images, including low to moderate noise levels and a range of realistic sampling kernels.

Algebraic curves, generalized moments, image sampling, signals with finite rate of innovation (FRI).

I Introduction

In today’s digital world, sampling is a key block of any signal acquisition device: the device senses and stores analog signals at certain points and uses the samples later for the representation of the analog signal (possibly after some post-processing). The main concern here is whether and how the collected samples provide a fair representation of the original signal. Hence, as a first step in the design of acquisition devices, we should develop suitable sampling and reconstruction techniques for the target class of signals.

The classical Shannon sampling theory and its variations present sampling strategies for bandlimited signals and more generally the class of signals living in a shift-invariant space [1, 2, 3, 4]. Still many crucial signals stay out of reach of this class. Among them are signals which can be described with a finite number of parameters, hence called signals with finite rate of innovation (FRI). In [5], a study of one-dimensional (1D) FRI signals was presented. This work then evolved to include more general FRI signals such as piecewise polynomials [5, 6], streams of Diracs [5, 6, 7] and piecewise sinusoids [8].

Extension of sampling schemes to images is an essential but challenging problem. Because of the sharp intensity transitions along edges, images are non-bandlimited. Also, the diverse geometry of the edges in typical images excludes them from the known shift-invariant spaces. Some preliminary efforts to generalize the FRI framework to images led to the sampling schemes with adequate sampling kernels for step-edge images and polygons [9, 10, 11]. In a recent work, an FRI-based sampling scheme is presented for images with more versatile edge geometries [12]. The curves in this model are zero level sets of a mask function that is a linear combination of a finite number of two-dimensional (2D) exponentials. The curve parameters are shown to satisfy an annihilation system of equations which could be solved directly, or more robustly as a minimization problem.

The curve model introduced in [12] is novel and further investigation is needed to reveal its descriptive power –i.e., the range of shape geometries and the number of free parameters required for generating a given shape in the range. On the other hand, a rich parametric model for 2D curves already exists in the literature: algebraic curves [13, 14, 15]. An algebraic curve is the zero level set of a bivariate polynomial of a finite degree. Algebraic curves can be decomposed into a finite number of smooth arcs. Nevertheless, they are dense, in the Hausdorff metric, among all smooth curves. Hence, every curve can be approximated by a sequence of algebraic curves arbitrarily closely [16]. This characteristic makes them an excellent candidate in modeling the general image boundaries.

We call a subset of the 2D plane with an algebraic boundary curve an algebraic domain and the restriction of it to the image plane an algebraic shape. According to a classical result [17], an algebraic domain of degree can be uniquely determined from its set of 2D moments of order less than or equal to . But as stated in [18], “there has been so far no constructive way of passing from the given moments to the unique algebraic domain, or equivalently to the defining polynomial”. In [18] and [19], the authors present an algorithm for the reconstruction of a subset of bounded algebraic domains –called quadrature domains– from their moments. However, moments are inherently very sensitive to noise and consequently, the suggested algorithm (as noted by the authors) suffers from severe numerical instabilities.

Moments have been used as the standard descriptors of 2D shapes in [20, 21, 22]. Also, there are some works on the exact calculation of moments of the shapes with parametric boundary curves in terms of the curve parameters, through nonlinear equations. Examples are [23] for polygonal shapes and [24] for shapes with wavelet and spline curves.

I-a Contribution

In this paper, we propose sampling and reconstruction techniques for algebraic shapes. We first derive a set of linear annihilation equations for the shape parameters with the coefficients being factors of 2D moments of the image. We prove that any solution of these equations will lead to a polynomial that vanishes on the boundaries of the original shape. By employing sampling kernels that reproduce polynomials like the well-known B-splines [25], we are able to calculate the shape moments from the samples.

Moments are inherently very sensitive to the noise and the reason is that noise in the image or the samples is boosted by polynomial factors before it contaminates the moments. To overcome this difficulty, we replace moments with some generalized moments that are still reproducible from the samples but do not amplify the noise. This is achieved by multiplying the monomials in the conventional moments with a function that is adjusted to the sampling kernel and decays at the image borders. The advantages are threefold: we get more stable moments that can be reproduced by a wider range of sampling kernels. Furthermore, we can extend our model to algebraic shapes with unbounded boundaries.

In any sampling problem, consistency of the reconstruction with noiseless samples is a crucial constraint [26, 27, 28]. It is also proved to be a strong tool for recovering binary images in the absence of a parametric model [29]. In this paper, we further improve the stability of our reconstruction by enforcing measurement (or sample) consistency to the recovered algebraic shape. This results in a reconstruction algorithm that is robust to moderate noise levels in the samples.

I-B Organization of the paper

The paper is organized as follows. In Section II, we first define the image model and study algebraic curves in details. Then, we explicitly define the sampling problem. We derive the annihilation equations for the shape parameters in Section III and present a perfect reconstruction algorithm for the noiseless scenario. In Section IV, we develop a stable reconstruction algorithm. For this purpose, we introduce the notion of generalized moments and present an algorithm for generating the adequate generalized moments corresponding to the given sampling kernel. Also, we prove that any solution of the annihilation equations formed from (generalized) moments generates the original shape boundaries. We present some experimental results with different curves in the noiseless and noisy scenarios in Section V and conclude in Section VI.

Ii Sampling of algebraic shapes

Ii-a Image model

Consider a bivariate polynomial of degree with real coefficients


The set of points defines an algebraic domain. The boundary of this domain, defined by the zero level set of , is an algebraic curve of degree ,

Let denote a closed domain in modeling the image plane. Without loss of generality, we take for some . We define an algebraic shape in as the binary image


where denotes the indicator function. This means that the edges of are contained in the algebraic curve .

An algebraic shape of degree is specified with parameters (the coefficients in (1)). In developing the annihilation equations of Section III, we assume that the algebraic shapes have closed boundaries. This restricts the polynomial degree to the even integers. We later remove this assumption by introducing generalized moments in Section IV.

Fig. 1: Algebraic domains of degree 4.

Typical examples of algebraic domains of degree 2 are circles and ellipses. Figure 1 displays two algebraic domains of degree 4. We see in this figure that an algebraic domain of degree 4 can have four disconnected components. The following remark asserts that this is an upper bound.

Remark 1 ([30]).

An algebraic domain of degree cannot have more than disconnected closed components.

This remark is a consequence of Bezout’s theorem [13]. We will also make use of this theorem in Section IV to prove our result.

Theorem 1 (Bezout).

Two algebraic curves of degree and that do not share a common component intersect in at most points.

Fig. 2: An algebraic shape of degree at least 4.

Bezout’s theorem also provides us a handy tool to roughly estimate the degree of an algebraic shape. Consider a shape image with boundary . should have a degree of at least if it intersects a line (a first-degree polynomial) at points or if it intersects an ellipse (a polynomial of degree 2) at points. This is illustrated with an example in Figure 2.

Algebraic curves have been studied and applied to data fitting in computer vision (e.g. [30, 31]). This rather long history of application has revealed that polynomials of modest degree (e.g. degree 4 with 15 parameters) have enough descriptive power to generate a diverse range of curve geometries. Hence, in the rest of this paper, we mostly consider . Nevertheless, all results remain valid for higher degree polynomials.

Ii-B Sampling

In a typical sampling setup (Figure 3), the image is first convolved with a 2D kernel and then sampled at a uniform grid to generate the samples

In a noisy setup, the noise vector will be added to the measurements after spacial sampling. The sampling kernel is determined by the physics of the sampling device but in most cases it can be considered as a separable kernel . In the first part of this paper, we consider separable kernels that can reproduce polynomials up to some degree. is a polynomial reproducing kernel of degree if there exist coefficients such that [32]

B-splines are well-known examples of polynomial reproducing kernels [25]. A zero order B-spline is defined as

A B-spline of order is obtained by convolving kernel

The B-spline kernel can reproduce monomials up to degree and the corresponding coefficients are obtained as

where is the dual of [25].

Fig. 3: In a typical sampling scenario, the image goes through convolution with a 2D kernel and spacial sampling to generate the measurements.
Fig. 4: (a) An algebraic shape. (b) Samples generated with the tensor product of B-spline kernels of order 6.

In any image sampling scenario, the question is whether and how we can reconstruct the original image from a finite number of samples (Figure 4). In the next sections, we present a technique for the reconstruction of the boundary curve and hence the algebraic shape from adequate noiseless or noisy samples .

Iii Reconstruction from Moments

For an exact reconstruction of an algebraic shape image, we should estimate its boundary –the algebraic curve – from the samples. In the sequel, we first derive some annihilating equations for the curve parameters based on the shape moments. Then, we use the existing FRI techniques [32] to calculate shape moments from the samples. The overall procedure is summarized in Algorithm 1.

Iii-a Annihilation equations

Consider a closed algebraic curve inside the domain and the corresponding shape image . We can rewrite in equation (2) as

where denotes the closure of the interior of . This equation explains that the partial derivatives and vanish everywhere in except possibly on , where they behave like the Dirac function. So, similar to the equation , we conclude that


inside .

We can multiply the above equations with for any and integrate over the domain to obtain the equations


By substituting from equation (1) in (5) and using integration by parts, we get


In the derivation of (7), we also used the fact that is a closed curve inside and hence, is zero at the domain borders.

Input: noiseless samples , degree of the algebraic shape, polynomial reproducing coefficients of the sampling kernel.
Output: boundary curve .

1:Calculate shape moments from samples for any , according to equation (III-A).
2:Form the annihilation equations (8) and (9) for any and put them into a linear system of the form .
3:Solve for the polynomial coefficients with the constraint .
4:Form the polynomial from the coefficients in according to (1).
5:Set equal to the zero level set of inside .
Algorithm 1 Algebraic shape reconstruction from noiseless samples

The integrals in equation (7) represent 2D moments of the image

Hence, we can rewrite equation (7) as


We can similarly modify equation (6) to derive the additional equation


For any pair of , formula (8) and (9) give us two linear annihilation equations for the coefficients , in terms of the image moments. We get enough equations to build a linear system of the form


and derive the curve parameters, if we consider all pairs . This implies that we require all image moments of degree up to , i.e., .

To avoid the trivial solution , we set the term corresponding to to 1. We recall that a scaling of the polynomial coefficients does not change its level sets. In Theorem 2, we prove that the zero level set of the polynomial formed by any solution of (10) contains . This specifically means that although the system of equations in (10) might have a null space with dimension larger than 1, any vector in this null space generates a polynomial that vanishes on the boundary of . Hence, we can recover the boundary curve and the algebraic shape from any solution of (10).

Fig. 5: The exponential growth rate of polynomial reproducing coefficients of the B-spline kernel .

Finally, it remains to retrieve moments from the samples. Suppose that the kernel can reproduce polynomials up to degree , with the corresponding coefficients . The 2D moments of the image can be calculated as


where and indicate the set of indices and such that is nonzero over .

Iii-B Stability

Algorithm 1 restores the exact algebraic curve when it has access to the noiseless samples. But it breaks down in the presence of noise. The reason is that the polynomial reproducing coefficients have the same growth rate as the polynomials, i.e., they grow like . (To illustrate this, we show the polynomial reproducing coefficients of a 1D 6th order B-spline kernel for in Figure 5.) This specially implies that in equation (III-A), the weight of samples that are away from the image center are considerably larger than the weight of the central samples. But for images in our model, samples at the image borders mostly contain noise. This transfers an amplified noise to the moments and results in severely degraded moments SNR. The noise boosting effect becomes more critical as the order of moments grow. This makes Algorithm 1 unstable even at a sample SNR as high as 100 dB.

We recall that in the related works of [33] and [11], only the first order moment are required as they focus on first degree polynomials (step edges). Hence, the aforementioned noise boosting effect is not an issue.

In the next section, we introduce some generalized moments that have slower growth rates and discard the noise at the image borders. Above all, they are still reproducible from the samples generated with a wider range of sampling kernels.

Iv Stable Recovery

The sampling scheme of Section III has some limitations: the reconstruction algorithm succeeds only in the absence of noise; the acceptable sampling kernels are limited to the ones that exactly reproduce polynomials; and the algebraic shapes should have closed boundary curves. In this section, we modify Algorithm 1 in three steps to resolve these limitations:

First and foremost, we introduce a fast decaying (or even compact-support) function in the integrands of equations (5) and (6) to reduce the growth rate of polynomials, especially at the borders of . This translates into the annihilation equations as replacing moments with some generalized moments. We prove in Theorem 2 that under noiseless samples, the resulting annihilation equations restore the exact boundary curve of any algebraic shape. Our proof is general and includes the case which leads to conventional moments. Next, we describe the requirements for to ensure stable generalized moments and we propose an optimization procedure for finding the best candidate that pairs with a given sampling kernel. Interestingly, the inclusion of allows for extension of the image model to algebraic shapes with open boundaries.

For our second step, we note that the image moments do not take full advantage of the available samples. For instance, the samples allow for prediction of the sign of the implicit polynomial on a subset of the sampling grid points and this prediction is fairly robust against noise. To further improve the reconstruction, we enforce sign consistency of the polynomial with the prediction of the available samples.

In our last step, we encourage full measurement consistency (not just sign) through bounded changes in the coefficients of the implicit polynomial.

Iv-a Annihilation equations with generalized moments

We developed the annihilation equations of Section III-A by multiplying equations (3) and (4) with . This caused the image moments to appear in the equations. To control the growth rate of the polynomials and hence the moments, we replace with for an appropriate function .

Definition 1.

For any bivariate function and integers , we call


a 2D generalized moment of , associated with .

Having separable sampling kernels, we also take to be separable of the form . Though, the results can be similarly extended to the non-separable kernels. In the following, we derive the new annihilation equations and discuss the requirements on afterwards.

We multiply equations (3) and (4) with and repeat similar steps as in section III-A to obtain


In Section III, we had to assume that is zero at the borders of the image plane in order to use integration by parts. Here, we assume that is either zero outside or decays so fast that the integral outside of this interval becomes negligible. This allows to take non-zero values at the borders of ; consequently, can represent an unbounded shape.

We can further simplify equations (13) and (14) and substitute the integrals with generalized moments to get the new annihilation equations


The above equations are valid for any . Note that restores the annihilation equations (8) and (9) when represents a closed shape. In Theorem 2 we state a unified result for recovery of algebraic shapes either from conventional annihilation equations or the generalizations in (15) and (16). The proof of this theorem is provided in the appendix.

Fig. 6: B-spline kernels and their associated ’s for reproducing stable generalized moments of order lass than or equal to 6. The indices () of the contributing kernels in equations (17) and (18) and the minimum number of required samples () are (a) and , (b) and , (c) and , respectively.
Theorem 2.

Let denote an algebraic shape of degree defined on without singular edges111We call an edge singular if the image color does not change on either of its sides; for instance the image associated with has a singular edge at points with equal coordinates.. Also, let and denote the generalized moments of (Definition 1) corresponding to a function for which vanishes outside and remains strictly positive inside . If satisfies the annihilation equations (15) and (16) for all , then, the zero level set of the polynomial

contains the boundaries (edges) of .

Remark 2.

Unique recovery of is not generally possible. Obviously, the zero level sets of and are the same, leading to the same algebraic shapes. However, there are less obvious examples that prevent unique recovery: the zero level sets of both and coincide with the unit circle, while the two bivariate polynomials have the same degree. The important point in Theorem 2 is that the curve is uniquely determined, but possibly with a different implicit polynomial.

Remark 3.

Theorem 2 requires that the coefficients of satisfy the annihilation equations (15) and (16) for every . This generates an over-determined system of the form with about 8 times more rows than columns. In our experiments, we have confirmed successful recovery of algebraic curves from the annihilation equations corresponding to (yielding an almost balanced system). Our proof technique, however, falls short of this stronger result.

Iv-A1 Optimal generalized moments

The primary reason of introducing to the equations is to control the growth rate of the monomials , especially at the image borders. Ideally, the function in (13) and (14) should be set such that and both vanish outside . The faster they decay near the borders of , the more stable will be the annihilation equations (15) and (16). However, the bottleneck in setting is the reproduction of moments from the samples. That is the functions and should be reproducible by the sampling kernel , i.e., we need coefficients and that satisfy


Here, represents values for which has an effective support in ; this ensures that and vanish outside .

For recovering an algebraic curve (domain) from samples using the generalized moment technique, we need to linearly combine the samples in correspondence to the coefficients and . In other words, we never require the function explicitly in practice. Consequently, instead of looking for the best function, we can search for coefficients and such that


To find such coefficients, we introduce the following objective function

where stands for the derivative of . Next, we solve the quadratic program


The equality constraint in the above minimization is to avoid the trivial zero solution and the inequalities guarantee that is non-negative. Although solving a quadratic program is computationally manageable, we have frequently observed that (20) is ill-conditioned222Essentially, the source is the same as the one causing instability in Algorithm 1 except there is no noise here: the error terms corresponding to different ’s in the objective function grow polynomially and this makes the problem ill-conditioned. in the sense that iterative methods are very slow in achieving the global solution, and usually terminate much earlier than desired. This shortcoming could be improved by using a sufficiently good initialization. Furthermore, any set of coefficients which result in a small cost according to the objective function could be used.

We recall that an implicit parameter in this problem is the size of the index set . This parameter also affects the modeling of and the minimum required sampling density for this sampling kernel. In fact, by increasing the index set the global cost in (20) can only reduce. Thus, the larger the the lower the cost. However, larger translates into more image samples, and consequently more complexity.

For the B-spline kernels, we found surprisingly good candidates that make the objective function almost zero. Figure 6 shows the kernels and their associated ’s that reproduce stable generalized moments of order 6 or less. This implies that we can form the annihilation equations and recover algebraic shapes of degree 4 even when the sampling kernel is the tensor product of 2nd order B-splines. The cost is a larger number of required samples.

Our final remark concerns using an asymmetric function in the form , when satisfying (17) and (18) is not possible simultaneously for a single . One can verify that the annihilation equations can still be obtained if

For finding and , (20) needs to be divided into two quadratic programs that accommodate and separately.

Iv-A2 Patch-based recovery

Equations (17) and (18) show that and consequently have compact support. This indicates that the generalized moments are computed from a finite window of the image samples –namely, of size , where amounts to the number of contributing kernels in . Having access to more samples, we can slide a window over the image samples, calculate 2D moments and form the annihilation equations for each window. This results in a linear system with more equations and improves the noise stability of the reconstruction.

Fig. 7: A compact support function facilitates the calculation of 2D generalized moments and the annihilation equations for different windows of the image samples. However, the coordinate shifts between different windows should be compensated before concatenating the equations in one system.

There is only one issue requiring further attention: in the annihilation equations of each window, the coordinates origin is taken at the window’s center (Figure 7). This means that the variables of each set of annihilation equations are the coefficients of the polynomial in those coordinates. Hence, we should compensate for the shifts in the coordinates before concatenating the equations of different windows. For this purpose, we choose the reference coordinates as the symmetry axes of the image plane. When the coordinates are shifted by , the polynomial in the original system shall be mapped to the polynomial

This reveals the mapping between the coefficients of , denoted by , and ’s as

for any . We can represent the above relations for all polynomial coefficients simultaneously as


where is an upper triangular square matrix with diagonal entries equal to 1. This allows us to relate the annihilation equations of a window centered at to the polynomial coefficients in the reference coordinate system through the equation

In a nutshell, we should multiply the annihilation equations of different windows with the corresponding matrix in equation (21) before concatenating them in a bigger system.

Iv-B Constraints on the sign of the polynomial

So far, we have built a system of equations in terms of the image parameters that is stable at numerical precision. In the presence of noise, the annihilation equations are only approximately singular. In this case, as a common practice, we consider the solution of the least squares minimization problem


The least squares denoising works well at low noise levels, especially when is a tall matrix. Nevertheless, since algebraic curves are dense among continuous curves, distortion in the image moments (originated from moderate noise levels in the samples) can lead to substantially different solutions.

Recently, the Cadzow’s denoising algorithm [34] has been used for denoising of the annihilation equations of 1D [7] and 2D [12] FRI signals. The common feature in these works that makes denoising successful is having annihilation equations with a Toeplitz structure. Our system of annihilation equations –although almost each element in has a few duplicates– is not Toeplitz and the Cadzow’s denoising algorithm does not help333In our implementation of Cadzow’s algorithm, we observed that it converges to a rank deficient matrix with the expected structure which stays very close to the noisy matrix ..

In our problem, the best reconstruction is an algebraic shape that is as consistent as possible with the image samples (i.e., up to the samples SNR). Theoretically, this can be achieved with a brute-force search over the space of image parameters. But this problem is nonconvex with many parameters and hence, computationally intractable. In the rest of this section, we exploit the local information provided by the samples to improve the reconstruction in the presence of noise.

Sample values represent the area of the intersection of the corresponding kernels with the interior of the shape in a weighted form. For example, indicates that everywhere in the support of 444We assume that and has a unit integral.. We further incorporate the samples in our reconstruction by interpreting them as the central points of the corresponding kernels lying inside or outside the shape. More precisely, if is above for an , we assume its center to be inside the shape, i.e., or equivalently . Also, we take or , if . Eventually, we constrain the solution of the least squares problem with the inferred signs:


where and encode respectively, the normal and sign-negated polynomial evaluation matrices at central locations of the sampling kernels; corresponds to locations with large sample values, while corresponds to locations with small sample values. The minimization problem (23) can be solved with quadratic programming algorithms.

Iv-C Measurement consistency

At moderate noise levels (sample SNRs around 25 dBs), the recovered curves from (23) are close enough to the original boundaries to let us approximate the function mapping the polynomial coefficients to the image samples with a 1st order Taylor expansion around the correct coefficients. We exploit this assumption to improve the measurement consistency of the reconstruction.

Let denote the mapping from the polynomial coefficients into the samples of the algebraic shape. For instance, if stands for the polynomial coefficients associated with an algebraic curve, represents the vector of noiseless image samples via the sampling kernel.

For a given set of noisy samples , let be the solution to the sign consistency technique in (23), which corresponds to . For moderate to low noise levels, we know that is a good approximation of . Thus, we use the linearization of around (1st order Taylor expansion) to write that

where is a matrix that relates the small input variations in to its output around the point . In practice, we find by numerically varying in all directions and observing the corresponding ’s. Finally, we improve our current estimate of by


In our algorithm, we apply few iterations of the above update rule. Each time we evaluate the associated vector and continue the iterations as long as this vector gets closer to .

V Experimental Results

We evaluate the performance of the proposed algorithm in different scenarios. We select bounded algebraic shapes for most of the experiments. For this purpose, we restrict the polynomial degree to even integers. But a randomly generated even degree polynomial very likely has unbounded level sets. A full characterization as well as a model for the generation of bivariate polynomials of degree 4 with bounded level sets was presented in [30] and [31]. We adopt this model to generate shapes for our experiments.

V-a Noiseless recovery

In the first experiment, we study reconstruction of algebraic shapes from noiseless samples. Recalling the results of the last section, we expect to recover the exact image by solving the least squares problem (22). Figure 8 displays perfect reconstruction of an algebraic shape of degree 4, when the sampling kernel is the tensor product of the 6th order B-splines.

Fig. 8: Exact reconstruction of algebraic shapes from noiseless samples. (a) An algebraic shape of degree 4. (b) Noiseless samples (size ), when the sampling kernel is . (c) Absolute difference between the original shape and the least squares solution.

V-B Recovery in the presence of noise

In this experiment, we aim at studying the effect of each step of the algorithm on the reconstructed image from noisy samples. For this purpose, we consider two distinct algebraic shapes of degree 4 with different levels of noise in their samples and we plot each stage of the reconstruction (Figures 9 and 10). The samples of both images are generated with the sampling kernel and the annihilation equations involve generalized moments corresponding to the function in Figure 6(a). We see that although the least squares solution might be offbeat in presence of noise, the constraints on the sign of the polynomial substantially restrain the solution and lead to satisfactory reconstructions at moderate signal to noise ratios (SNRs).

Fig. 9: Reconstruction from noisy samples. (a) Noisy samples of the shape in Figure 8(a) with size and SNR = 17 dB. (b) Absolute error of the least squares solution (PSNR = 13.7 dB). (c) Ablsoute error of the quadratic programming (equation (23)) recontsruction (PSNR = 20.4 dB). (d) Absolute error of the output of the consistency improvement algorithm (PSNR = 21.3 dB). SNR between the samples of the final reconstruction and the noisy samples (a) is 15.4 dB.
Fig. 10: Reconstruction from noisy samples. (a) Original image. (b) Noisy samples of size with SNR = 22 dB. (c) Absolute error of the least squares solution. (d) Absolute error of the quadratic programming solution (PSNR = 21.0 dB). (e) Absolute error of the output of the consistency improvement algorithm (PSNR = 22.8 dB). SNR between the samples of the final reconstruction and the noisy samples (b) is 20.9 dB.

V-C Sampling kernel sensitivity

Earlier, we mentioned that a consequence of replacing conventional moments with generalized moments is relaxing the restrictive polynomial-reproducing requirement on the sampling kernel. Specifically, we worked out the reproducing coefficients for the B-spline kernels of order 2, 4, and 6 that generate stable generalized moments of order less than or equal to 6 (see Figure 6). This, for example, allows us to recover algebraic shapes of degree 4 from samples generated with the sampling kernel . In this experiment, we study the sensitivity of the reconstruction to the choice of the kernel. Figure 11 displays the absolute difference between an image and its reconstructions from samples generated with different sampling kernels and similar signal-to-noise-ratios. The results are comparable irrespective of the choice of the sampling kernel (note the expected difference in the sample sizes that calls for different noise realizations for the three samples).

Fig. 11: Sensitivity of the reconstruction to the choice of the sampling kernel. (a) Algberaic shape of degree 4. (b),(c),(d) Noisy samples (SNR = 27 dB) of size , and , generated with B-spline kernels of degree 2, 4, and 6, respectively. (e),(f),(g) Absolute error of the reconstructions from samples in (b),(c) and (d) with reconstruction PSNRs 22.2 dB, 23.2 dB and 21.3 dB, respectively.

V-D Unbounded algebraic shapes

Introducing generalized moments to the annihilation equations facilitated sampling and reconstruction of algebraic shapes with open boundaries (also referred as unbounded algebraic shapes). This additionally allows the reconstruction to enjoy oversampling by forming annihilation equations for each sample window, without caring about the image content of the window. Figure 12 shows the reconstruction of an unbounded image from its noisy samples, where the sampling kernel is the tensor product of 2nd order B-splines. The peak signal-to-noise ratio (PSNR) of the reconstructed image is 20.1 dB and SNR between its samples and the original noisy samples (sample consistency) is equal to 23.1 dB. These numbers clearly indicate the success of our proposed algorithm for reconstructing unbounded shapes.

Fig. 12: (a) Unbounded algebraic shape of degree 4. (b) Noisy samples of size with SNR = 25 dB. (c) Absolute reconstruction error.

V-E Overfitting

In the last two experiments, we address uncertainties in the image model. First, we study the situation where we overestimate degree of an algebraic shape. Recalling Theorem 2 of the previous section, we expect the recovered polynomial from the annihilation equations to vanish on the boundaries of the original shape in the noiseless scenario. Figure 13 displays the results when we approximate an ellipse with algebraic shapes of degree 4 from its noiseless and noisy samples. Figures 13(c) and 13(f) show the least squares solutions for noiseless and noisy samples, respectively. Both figures indicate that the boundaries of the recovered images contain the boundary of the original ellipse. Equivalently, the recovered polynomials are factors of the original polynomial of degree 2. The extra factors are resolved in the next steps of the algorithm, resulting in exact reconstructions in both scenarios.

Fig. 13: Approximation of an ellipse with algebraic shapes of degree 4. (a) Original ellipse. (b) Noiseless samples of size , generated with B-splines of degree 2. (c) Least squares solution for noiseless samples. (d) Absolute error of the final reconstruction of the algorithm. (e) Noisy samples with SNR = 22 dB. (f) Least squares solution for noisy samples. (g) Absolute error of the final reconstruction from noisy samples.

V-F Algebraic shape approximation

Another type of uncertainty in the image model happens when the image boundary is not an algebraic curve. Regarding the descriptive power of algebraic curves, we still expect to find a good approximation of the image. To investigate this, we generated a shape with a Bézier curve boundary with four control points and generated its -samples with 2nd order B-spline sampling kernels. Then, we obtained the approximate algebraic shape from the noiseless samples. The original image and the absolute error of its algebraic approximation are depicted in Figure 14. We observe that the reconstructed curve is a rather accurate descriptor of the original Bézier curve.

Fig. 14: Approximation of non-annihilable curves with algebraic curves. (a) A shape with a Bézier curve boundery. (b) noiseless samples. (c) The absolute error between the original shape and its approximation with an algebraic shape of degree 4. The reconstruction PSNR is 19.8 dB.

Vi Conclusion

Designing sampling schemes for images with arbitrary edge geometries is still a challenging research problem. In this paper, we proposed a sampling and reconstruction algorithm for binary images with boundary curves that are zeros of an implicit bivariate polynomial. We developed a set of linear annihilation equations from the image samples and proved that every solution of the equations restores the image boundaries, in the noiseless scenario. The primary equations involve 2D moments of the image. To make the reconstruction robust against noise, we replaced conventional moments with generalized moments associated with a compact-support 2D function that is paired with the given sampling kernel. This leads to a reconstruction algorithm from more realistic samples and extends the model to images with open boundaries.

The image model we considered in this paper is very rich and may be used for the approximation of general shapes from their samples. Also, the idea of replacing conventional moments with generalized moments might find applications in other image processing tasks which use moments as the image descriptors.


The authors would like to thank P. Thévenaz for his helpful hints and suggestions. Also, they are thankful to L. Baboulaz for sharing MATLAB codes and P. L. Dragotti for the valuable discussion in the beginning of this work.

[Proof of Theorem 2] We prove by contradiction. Assume the zero level set of does not fully include ; thus, can be factorized as

where is coprime with and , and has a non-trivial zero level set . Meanwhile, the zero level set of is included in that of . Roughly speaking, and stand for parts of that are excluded and included in the zero level set of , respectively. Further, let be a polynomial with minimum degree such that . If is irreducible, we shall have , otherwise, might be different from . In either case, we have .

The validity of annihilation equations (15) and (16) imply

for all . By linearly combining these equalities, we conclude that


holds for any polynomial of degree no higher than . From this point on, we set as


Because , this choice of fulfills the degree constraint.

Let be such that the line intersects . According to Bezout’s theorem, the number of intersections shall be limited to . We assume the intersections are at and conclude that


where is the Dirac’s delta function and are sign values; () if is positive (negative) at and negative (positive) at for small enough . Hence,

This shows that the value of at is either or has the opposite sign as . This implies that


where equality happens only if . By taking advantage of (27), we can rewrite the inner integral in (25) as


Thus, for (25) to hold, needs to vanish at all points on , and in particular, at points on . As and are coprime, can vanish only on a finite number of points on (Bezout’s theorem) . This forces the zero level set of to include (inclusion of except finitely many points implies inclusion of the whole ).

For any , because of we have that


Again, since and are coprime, can happen only for a finite number of points . Therefore, should hold for all ; i.e., the zero level set of includes the zero level set of . This, however, contradicts our initial assumption that is a polynomial with minimum degree that satisfies this property.


  • [1] C. E.Shannon, “A mathematical theory of communication,” Bell. Syst. Tech. J., vol. 27, pp. 379–423, 1948.
  • [2] A. J. Jerri, “The Shannon sampling theorem: Its various extensions and applications: A tutorial review,” Proc. IEEE, vol. 65, no. 11, pp. 1565–1596, Nov. 1977.
  • [3] M. Unser, “Sampling-50 years after Shannon,” Proc. IEEE, vol. 88, no. 4, pp. 569–587, Apr. 2000.
  • [4] M. Vetterli, J. Kovacević, and V. K. Goyal, Foundations of signal processing.   Cambridge University Press, 2014.
  • [5] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE Trans. Signal Process., vol. 50, no. 6, pp. 1417–1428, Jun. 2002.
  • [6] I. Maravić and M. Vetterli, “Sampling and reconstruction of signals with finite rate of innovation in the presence of noise,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 2788–2805, Aug. 2005.
  • [7] T. Blu, P. L. Dragotti, M. Vetterli, P. Marziliano, and L. Coulot, “Sparse sampling of signal innovations,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 31–40, Mar. 2008.
  • [8] J. Berent, P. L. Deragotti, and T. Blu, “Sampling piecewise sinusoidal signals with finite rate of innovation methods,” IEEE Trans. Signal Process., vol. 58, no. 2, pp. 613–625, Feb. 2010.
  • [9] I. Maravić and M. Vetterli, “Exact sampling results for some classes of parametric nonbandlimitted 2-d signals,” IEEE Trans. Signal Process., vol. 52, no. 1, pp. 175–189, Jan. 2004.
  • [10] P. Shukla and P. L. Dragotti, “Sampling schemes for multidimensional signals with finite rate of innovation,” IEEE Trans. Signal Process., vol. 55, no. 7, pp. 3670–3686, Jul. 2007.
  • [11] C. Chen, P. Marziliano, and A. C. Kot, “2D finite rate of innovation reconstruction method for step edge and polygon signals in the presence of noise,” IEEE Trans. Signal Process., vol. 60, no. 6, pp. 2851–2859, Jun. 2012.
  • [12] H. Pan, T. Blu, and P. L. Dragotti, “Sampling curves with finite rate of innovation,” IEEE Trans. Signal Process., vol. 62, no. 2, pp. 458–471, Jan. 2014.
  • [13] R. J. Walker, Algebraic curves.   Princeton, New Jersey, 1950.
  • [14] W. Fulton, Algebraic curves: an introduction to algebraic geometry.   Addison-Wesley, 1989.
  • [15] E. Arbarello, M. Cornalba, and P. Griffiths, Geometry of algebraic curves: volume II with a contribution by Joseph Daniel Harris.   Springer Science & Business Media, 2011, vol. 268.
  • [16] B. Gustafsson, “Quadrature identities and the Schottky double,” Acta. Appl. Math., vol. 1, no. 209–240, 1983.
  • [17] S. Karlin and W. J. Studden, Tchebycheff systems, with applications in analysis and statistics.   New Yourk: Interscinece, 1966.
  • [18] P. Milanfar, M. Putinar, J. Varah, B. Gustafsson, and G. Golub, “Shape reconstruction from moments: theory, algorithms, and applications,” in International Symposium on Optical Science and Technology.   International Society for Optics and Photonics, 2000, pp. 406–416.
  • [19] B. Gustafsson, C. He, P. Milanfar, and M. Putinar, “Reconstructing planar domains from their moments,” Inverse Probl., vol. 16, no. 4, pp. 1053–1070, 2000.
  • [20] S. Rad, K. C. Smith, and B. Benhabib, “Application of moments and Fourier descriptors to the accurate estimation of elliptical shape parameters,” in Proc. 1991 IEEE Int. Conf. Acoust. Speech, Signal Process. (ICASSP), vol. 4, 1991, pp. 2465–2468.
  • [21] R. Desai, R. Cheng, and H. D. Cheng, “Pattern recognition by local radial moments,” in Proc. 12th IAPR Int. Conf. Pattern Recog. (ICPR), 1994.
  • [22] K. Tsirikolias and B. G. Mertzios, “Statistical pattern recognition using efficient two-dimensional moments with applications to character recognition,” Pattern Recog., vol. 26, no. 6, pp. 877–882, Jun. 1993.
  • [23] M. Singer, “A general approach to moment calculation for polygons and line segments,” Pattern Recog., vol. 26, no. 7, pp. 1019–1028, Jul. 1993.
  • [24] M. Jacob, T. Blu, and M. Unser, “An exact method for computing the area moments of wavelet and spline curves,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 23, no. 6, pp. 633–642, Jun. 2001.
  • [25] M. Unser, “Splines: a perfect fit for signal and image processing,” IEEE Signla Process. Mag., vol. 16, no. 6, pp. 22–38, Nov. 1999.
  • [26] N. T. Thao and M. Vetterli, “Deterministic analysis of oversampled A/D conversion and decoding improvement based on consistent estimates,” IEEE Trans. Signal Process., vol. 42, no. 3, pp. 519–531, Mar. 1994.
  • [27] M. Unser and A. Aldroubi, “A general sampling theory for nonideal acquisition devices,” IEEE Trans. Signal Process., vol. 42, no. 11, pp. 2915–2925, Nov. 1994.
  • [28] M. Unser and J. Zerubia, “A generalized sampling theory without band-limiting constraints,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 45, no. 8, pp. 959–969, Aug. 1998.
  • [29] M. Fatemi, A. Amini, L. Baboulaz, and M. Vetterli, “Shapes from pixels,” arXiv preprint arXiv:1508.05789, 2015.
  • [30] D. Keren, D. Cooper, and J. Subrahmonia, “Describing complicated objects by implicit polynomials,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 16, no. 1, pp. 38–53, Jan. 1994.
  • [31] G. Taubin, F. Cukierman, S. Sullivan, J. Ponce, and D. J. Krigman, “Parameterized families of polynomials for bounded algebraic curve and surface fitting,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 16, no. 3, pp. 287–303, Mar. 1994.
  • [32] P. L. Dragotti, M. Vetterli, and T. Blu, “Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix,” IEEE Trans. Signal Process., vol. 55, no. 5, pp. 1741–1755, May 2007.
  • [33] L. Baboulaz and P. Dragotti, “Exact feature extraction using finite rate of innovation principles with an application to image super-resolution,” IEEE Trans. Image Process., vol. 18, no. 2, pp. 281–298, Jan. 2009.
  • [34] J. A. Cadzow, “Signal enhancement–a composite property mapping algorithm,” IEEE Trans. Acoust. Speech, Signal Process., vol. 36, no. 1, pp. 49–62, Jan. 1988.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description