Stable recovery of planar regions with algebraic boundaries in Bernstein form

We present a new method for the stable reconstruction of a class of binary images from a small number of measurements. The images we consider are characteristic functions of algebraic domains, that is, domains defined as zero loci of bivariate polynomials, and we assume to know only a finite set of uniform samples for each image. The solution to such a problem can be set up in terms of linear equations associated to a set of image moments. However, the sensitivity of the moments to noise makes the numerical solution highly unstable. To derive a robust image recovery algorithm, we represent algebraic polynomials and the corresponding image moments in terms of bivariate Bernstein polynomials and apply polynomial-generating, refinable sampling kernels. This approach is robust to noise, computationally fast and simple to implement. We illustrate the performance of our reconstruction algorithm from noisy samples through extensive numerical experiments. Our code is released open source and freely available.


Introduction
In this work, we present an improved method for the sampling and recovery of a class of binary images of the form I = χ D , where D ⊂ R 2 is a bounded open region whose boundary ∂D is an algebraic curve of degree n, i.e., the zero locus of a bivariate polynomial p of degree n: We refer to such region D as an algebraic domain or algebraic shape and, without loss of generality, we assume that it is contained inside the rectangular region Ω = [0, L 1 ] × [0, L 2 ] ⊂ R 2 . By classical results [12], an algebraic curve of degree n as above can be uniquely determined from its set of two-dimensional moments of order less than or equal to n. An algorithmic approach for the reconstruction of bounded algebraic domains from their moments was first presented in [10,19] but this approach is very sensitive to noise. Even though it was shown that one can improve the stability of the reconstruction by increasing the number of moments [15], the problem of recovering an algebraic domain remains unstable in the sense that small errors in the computation of the image moments may have a significant impact on the recovery algorithm.
In this work, we adopt the setting recently proposed in [9] where it is assumed that we have access only to a discrete set of uniform samples of the binary image I = χ D , that is, input data are of the form where φ : R 2 → R is an appropriate sampling kernel and T ≥ 1 is a parameter. The problem we consider is whether we can recover the exact image I and, hence, the corresponding algebraic curve ∂D from an adequate set of noiseless or noisy samples {d k }, given the sampling kernel φ and the sampling rate T , where T may be large, that is, the image I may be heavily down-sampled. As in [9], we assume that φ is polynomial generating up to a certain degree m, that is, there exist coefficients c (α) k such that Fatemi et al. [9] show that, under the assumption stated above, one can express the image moments (2) as appropriate linear combinations of the image samples (3), where the coefficients of the linear combinations depend on the sampling kernel φ. Hence, they derive a system of linear equations for the recovery of the algebraic shape of the form M a = 0, where a is the vector of the unknown polynomial coefficients {a i,j } in (1) and the matrix M consists of the computed moments (2). It turns out that the direct numerical solution of (5) recovers the polynomial coefficients if the image moments are computed from noiseless image samples. However, if the image samples are corrupted by even a small additive noise, then the numerical reconstruction fails in general.
To remedy the instability of the recovery, Fatemi et al. [9] introduced a modified formulation based on "generalized moments" leading to a constrained optimization problem that is solved using an iterative regularized reconstruction algorithm. The instability of the solution of the system (5) can be explained by observing that the process of converting image samples into image moments can be very sensitive to noise. We contend that this sensitivity is highly dependent on the basis selected for the representation of the algebraic curve and so the impact of the noise can be reduced by choosing alternative polynomial representations. Hence, to derive a more robust method for the recovery of an algebraic curve from the corresponding image samples, we introduce a novel approach that represents an algebraic curve in terms of non-separable bivariate Bernstein polynomials. It is already known that representations in terms of Bernstein polynomials can provide enhanced stability in problems from numerical analysis (cf. [8,22]) and our results are consistent with this general observations. Using our representation of algebraic domains in terms of Bernstein polynomials, we derive a new formulation of the image moment equation (5) that we can solve directly to recover an algebraic domain from noisy image samples. We show that the numerical reconstruction based on our algorithm is robust to noise and, while its reconstruction performance is comparable with the best regularized algorithm in [9], it is computationally much faster and simpler to implement. This is particularly true in the situation where we use refinable sampling kernels which are polynomial reproducing due to the costless computation of the expansion coefficients in (4). In addition, our method is more flexible as it applies to algebraic shapes of any degree and to image samples with different noise levels without the need for noise-specific parameter tuning that is required by a regularization approach.
To provide a broader context to this work, we recall that the problem of accurately recovering image boundaries or edges from image samples has gained renewed interest in recent years with the study of signals with finite rate of innovation (FRI). This area of investigation is concerned with signals that are not band-limited, hence do not satisfy the assumptions of Shannon sampling theory; however, they can be described with a finite number of parameters so that they can still be recovered from their samples using appropriate alternative strategies [2,17,27]. In an effort to apply the FRI framework to images, a number sampling schemes with different sampling kernels were proposed to recover special classes of images with edges [3,20,23,29]. While prior results were mostly focused on images with polygonal shapes or edges associated with finite sums of complex exponentials, the more recent work in [9] was the first to consider a rather general class of binary images associated with algebraic domains. We also recall that the application of image moments to signal processing and pattern recognition has a long history with a number successful applications to problems including shape analysis for aircraft identification, scene matching, character recognition, landmark detection and image retrieval [11,21,24,26,30].
The rest of the paper is organized as follows. In Section 2, we recall some basic properties of Bernstein polynomials and formulate the problem of the recovery of an algebraic curve in Bernstein form from the solution a set of moment equations. In Section 3 we formulate the computation of the Bernstein moments in terms of image samples (3). Finally, in Section 4 we present numerical experiments to illustrate the performance of our algorithm in reconstructing binary images of algebraic domains from noisy, including low to moderate noise levels, different down-sampling rates and a range of sampling kernels. Numerical results show that the reconstruction performance of our approach is very competitive with respect to the regularized reconstruction in [9] at a much lower computational cost.

Reconstruction from Bernstein moments
In this section, we derive a Bernstein representation of a bounded open algebraic where p is a bivariate polynomial of degree n ≥ 1. For such a domain, the boundary ∂D is the algebraic curve It is known that there are many polynomials producing the same boundary. In the following, we will show that it is possible to establish conditions so that the moments of D can be used to derive a unique polynomial p associated with its boundary We assume that D is contained inside a rectangular region Ω = [0, L 1 ] × [0, L 2 ]. Since non-separable Bernstein polynomials (see [8,14]) are defined on triangular domains, rather than considering Ω as the union of two triangles sharing a common edge (and thus considering a piecewise Bernstein representation of p), we adopt here the simpler and less expensive solution of embedding Ω into a triangle T L with vertices (0, 0), (L, 0), (0, L), where L = L 1 + L 2 , as illustrated in Fig. 1.
We then can write where The n+2 2 free parameters for the Bernstein representation of p in Ω ⊂ T L can be obtained by considering the binary image I = χ D . To this purpose, we follow the arguments from [15] in order to derive conditions for the boundary ∂D represented in terms of Bernstein polynomials. We start by recalling the Stokes' theorem according to the generalization due to Whitney [28].
Theorem 1 Let D ⊂ R 2 be bounded and open, with boundary ∂D that is smooth up to a set of measure zero in R. Let n(x) be the outward pointing normal to D at x ∈ ∂D. Then, given a vector field X on R 2 and a differentiable function f , we have We next take f (x) = x α p(x), where α ∈ N 2 and p is a polynomial of degree n ≥ 1 vanishing on ∂D. Letting X = x = (x 1 , x 2 ) (vector field in R 2 ), Stokes' theorem gives the following: We now consider that p is given as in (6) and make use of some straightforward properties of the Bernstein polynomials (for details, see [14]) recalled in the following propositions. The first property is related to the multiplication of a Bernstein polynomials with a monomial factor. The second result concerns the partial derivatives of Bernstein polynomials. Proposition 1 For a fixed α ∈ N 2 , |α| = α 1 + α 2 , the following equality holds:

Proposition 2
The partial derivatives of the polynomial B n i,j with respect to the variables x 1 , x 2 are respectively given by with the convention that B n i,j (x) = 0 whenever i < 0, i > n, j < 0, j > n.
From the expressions above in particular we derive which, in virtue of (9), takes the form where we have introduced the Bernstein moments (B-moments) to be computed: Defining the matrix elements fixing l ∈ N 0 , from (11) and (13) we derive the homogeneous system of linear equations Setting b = (1, b 1,0 , b 0,1 , . . . b n,0 ) T and s 2 (l) = l+2 2 , the system (14) can then be written as . By taking a particular value of l, that is, by specializing the variation of α ∈ N 2 l , the non trivial solution to (11) yields the searched polynomial coefficients.
Existence and uniqueness of such polynomial as a solution of (14) is a direct consequence of Theorem 2.2 in [15] combined with a change of basis argument. We state below a version of the same theorem specialized to our setting.

Theorem 2 Let D ⊂ R 2 be a bounded open set with real algebraic boundary.
Assume that D = intD, the boundary ∂D has total degree n and the point x = 0 does not belong to the zero set of the ideal I (∂D). Let G n l ∈ R s 2 (l)×s 2 (n) be the matrix of coefficients given by (14) and associated with the moments of D. Then the system of linear equations Remark 1 Under the assumptions of Theorem 2, there is a unique set of polynomial coefficients solving the moment equation G n l b = 0, where uniqueness is up to normalization of the leading coefficient of the polynomial. The hypotheses of Theorem 2 exclude images whose boundaries possess singular points such as boundaries represented by the zeros of the polynomial p( In such situations, neither existence nor uniqueness of a polynomial as a solution of (11) can be assured (see the proof of Theorem 2.2 in [15]).
We conclude the section by noticing that, re-arraging the terms in (11) and using the equalities Therefore, a solution of (11) can be found by searching for a common solution of two systems We remark that [9] uses a less formal approach without invoking Stokes' Theorem to derive a system of linear equations similar to (15) in the simpler setting of a standard power basis.

Computation of the Bernstein moments
In order to set up the equation (11), we need a tool for the computation of the Bmoments (12) from the image samples (3). First of all we observe that, from our assumption on I , we have We assume that the kernel φ possesses the polynomial generation property up to degree n + |α|, so that there exists a set of coefficients {C Then the B-moments simplify as and, in virtue of (3), for T = 1, they reduce to Their computation then just requires the Bernstein polynomial expansion coefficients C i,j, k . A strategy for determining them is reported in the following subsection.

Computation of the expansion coefficients
It is clear that the computation of the Bernstein moments strongly depends on the selected kernel φ. An appropriate choice for φ is a refinable function, satisfying the refinement equation for some mask a = {a k } k∈Z 2 . Typically φ can be chosen as generating function of a multiresolution analysis (if it has stable or orthonormal integer translates) and can be seen as the limit function of a convergent subdivision scheme. Typical examples include B-splines and Daubechies scaling functions (see for example [4] for comprehensive references). We recall that a bivariate (binary) subdivision scheme consists of the repeated application of a linear subdivision operator S a : (Z 2 ) → (Z 2 ), associated with a bivariate sequence mask a = {a k } k∈Z 2 transforming, at each step k, a sequence p (k) of points in R 2 into a refined sequence of points in  The scheme is called convergent if for any initial sequence p (0) there exists a continuous function f p (0) (f p [0] ≡ 0 for at least one initial sequence The limit is often denoted by S ∞ (p (0) ).
The connection with refinable functions in (17) is that the limit obtained starting with the delta-sequence δ = {δ 0,α } i∈Z 2 , φ a := S ∞ (δ), usually called the basic limit function of the scheme, is indeed refinable with mask given by the subdivision coefficients (see [7] for details about subdivision schemes or the recent paper [6]).
If we choose the kernel φ to be the basic limit function of a subdivision scheme generating polynomials up to degree m, for example tensor product B-splines or Box-splines, the determination of the coefficients of the expansion of any polynomial p ∈ Π m in terms of translates of φ can be efficiently based on subdivision. Note that, in the lucky case the coefficients (of the expansion of any polynomial p ∈ Π m in terms of translates of φ) are simply the values of the polynomial at integers (possibly shifted), the kernel φ is said to be reproducing polynomials up to degree m. In what follows, we illustrate the determination of the coefficients in the fist situation, starting from the idea proposed in [16].
We make use of the symbol p referring both to the polynomial and the sequence of integer samples of the polynomial p = {p(α), α ∈ Z 2 }. The difference will be understood from the context.
For p ∈ Π m we consider its Taylor expansion and recall that since a polynomial is uniquely defined by its values on the integer grid and in consideration of the generation properties of φ, we have Therefore, using (18) and setting M α = ∈Z 2 φ(− ) α for α ∈ Z 2 , we obtain Now, if we search for the coefficients C p α , α ∈ Z 2 , of the expansion of any polynomial p ∈ Π m in terms of translates of φ, i.e., by identifying the right hand side of (19) as the action of the linear operator S ∞ (which, as denoted above, is associated to the repeated application of the subdivision operator), on Π m , we can search for the "inverse" of S ∞ on Π m , i.e., the operator Q satisfying Indeed, if p ∈ Π m we have so that the expansion coefficients are found as C p α = (Qp)(α), α ∈ Z 2 . To identify Q, we require that it can be written as: where the coefficients A α , |α| ≤ m, are solution of the system of linear equations obtained by imposing conditions (21), as shown below. That is, for p ∈ Π m we have or, equivalently, after setting α + β = γ and denoting I m = { ∈ Z 2 : | | ≤ m} and I γ m = { ∈ I m : ≤ γ }), Next, observing that if |γ | ≥ m then D γ p = 0 (since p ∈ Π m ) and writing I 2m as the union of I m and its complement, we have This is equivalent to the conditions The resulting system of linear equations N a = e 1 where a = (A (0,0) , A (1,0) ...., A (n,0) , A (0,1) , · · · , A (n−1,1) , · · · , A (0,n) ) T ( !) 2 for = 0, · · · ,m−1, wherem = (m+1)(m+2) 2 . Therefore, multiplication by the diagonal matrix yields to the equivalent system of linear equations (PN )a = e 1 with diagonal entries equal to 1. As a consequence, there exists a norm · such that a ≤ 1.
In the case p = B i,j , the coefficients of the generation of Bernstein polynomials will then be given by where the operator Q involves the derivatives D α B i,j . It is important to notice that, since these coefficients can be expressed in terms of Bernstein polynomials of lower degrees (see Proposition 2), then the amplitude of such coefficients does not depend on the size of the domain. This is a very significant difference with the respect to the situation where a power basis is used (see Section B in [9]).

Explicit computation of the coefficients through the univariate case
An important advantage of embedding the image domain Ω in the larger triangle T L , is the possibility of writing the bivariate Bernstein basis in a separable way as the product of two univariate Bernstein polynomials, as stated in the following proposition. This result, together with the assumption of separability on the kernel φ, gives us the possibility of finding the coefficients of the generation of the Bernstein polynomials required for computing (16) using just univariate techniques.

Proposition 3 The following relation holds:
So, let us impose the following assumptions on the kernel φ: 2. ϕ is a refinable function, satisfying the refinement equation 3. ϕ either generates or reproduces polynomials up to the degree m, that is there exist sequences c := {c k } k∈Z such that for every univariate polynomial p of degree ranging from 0 to m.
The number of coefficients involved in (25) depends on the support of the function ϕ and on the interval where t varies. Let supp ϕ = [0, N] and t ∈ [a, b], with a, b ∈ Z. Then c k is different from zero only for a + 1 − N ≤ k ≤ b − 1 so that the total number of coefficients is b − a + N − 1. Clearly, such ϕ can also reproduce or generate the Bernstein polynomials, up to the degree m. Using (22), we then obtain the following useful simplification of (16) on the triangle T L . Observing that it follows that we can express the coefficients in (16) as In conclusion, everything reduces to the computation of the coefficients in the univariate case. Therefore, we continue by specializing the algorithm described in the previous section to this situation. We have where M j are the discrete moments of ϕ: It is important to mention that such moments can be computed even if we do not have the analytic expression of the kernel ϕ. In fact, since ϕ is refinable, the values at the integers can be found by solving an eigenvector problem obtained by evaluating (24) at the integers. Setting i + j = k, we write which is equivalent to imposing the condition The resulting system of linear equations N a = e 1 is upper triangular with diagonal entries given by N , = M 0 ! for = 0, · · · , m and can be easily solved by backward substitution.
From the coefficients A i , i = 0, · · · , m (that are dependent on the refinable function ϕ) we obtain the coefficients of the expansion of any univariate polynomial p of degree m in the generic interval [a, b] in terms of the refinable function ϕ; that is showing that the coefficients c m k depend on the derivatives of the polynomial p evaluated at the integers. In the case of a Bernstein polynomial B m j on the interval [a, b], the ith derivative is given by: A comparison with the derivatives of the monomials allows us to conclude that a representation in terms of Bernstein polynomials rather than a power basis is much more stable due to the limited growth of the coefficients involved in (16).
Numerical evidence of such aspect is given in Fig. 2 where the Bernstein vs. the monomial coefficients are plotted, with respect to the following kernels: Daubechies 7, B-spline of order 8 (both with polynomial generation order 8) and the "dual" pseudospline (see [5]) possessing the polynomial reproduction property of the same order. The interval here is fixed to [−10, 10], even though the displayed coefficients lay on a different range, which takes into account also the support of the kernels as explained above.
Finally, we can summarize the steps of our image recovery approach as follows. As above, we assume that we are given image samples (d k ) computed according to (3) using a bivariate kernel φ associated with a refinable function ϕ as in (23)

Numerical validation
To validate our improved method for the recovery of a digital image I of an algebraic shape from samples {d k } of the form (3), we ran extensive numerical experiments where we considered different noise levels, sampling rates and sampling kernels. To benchmark our results, we compared the performance of our Bernstein-moments reconstruction (BMR) method against the conventional moment reconstruction (CMR) approach that uses a power basis in the formulation of the moment equations and the regularized moment reconstruction (RMR) algorithm that is proposed in [9].
In accord with our theoretical framework presented in Section 2, we selected bounded algebraic shapes for our numerical experiments. In particular, similar to [9], we restricted our examples to stably bounded polynomials, a subclass of bivariate polynomials with bounded level sets that are characterized in [13] and [25]. Even though our recovery algorithm applies to shapes associated with algebraic curves of any degree, for simplicity we considered algebraic curves of degree 4 for most of our examples.
To set up our numerical experiments, we partitioned the planar region [0, 20] × [0, 20] uniformly. Next, we randomly generated stably bounded polynomials p : R 2 → R such that their pre-image S p ⊂ R 2 of (−∞, 0] is a subset of [0, 20]×[0, 20] and, correspondingly, obtained images I p = χ S p of size 512 × 512. Here, there is no loss of generality about the parameters 20 and 512, and they can be changed without affecting the algorithm. To generate the samples {d k } we used B-spline sampling kernels as a default setting; we discuss the impact of changing the sampling kernel below.
Our numerical code was implemented using MATLAB Release 2019a [18] and is designed to recover the Bernstein polynomial coefficients of an algebraic curve from a set of image samples following the steps summarized in our Algorithm 1. Our code uses the MATLAB function pinv to solve the system of linear equations in Algorithm 1 as a Moore-Penrose pseudo-inverse and is available in GitHub. 1 Our numerical experiments were run on a 64-bit Fedora 30 workstation with a 4×Intel®Core™ i3-4130 CPU @ 3.40 GHz and 15.5 GB of RAM.

Image recovery experiments
In the absence of noise, our BMR algorithm recovers images of algebraic shapes from their samples very accurately and there is no significant difference in using a power basis or Bernstein polynomials for the formulation of the moment equations (15). However, as observed above, the CMR algorithm that uses a power basis to formulate the moment equations is very sensitive to noise and its performance degrades very sharply even for very low level of noise as we show below. As remarked in [9] and in our discussion above, this instability is due to the behavior of the polynomial reproducing coefficients c i k that exhibit the same growth rate as the corresponding polynomial basis, i.e., they grow like |k| i when the polynomial is represented as a power basis. To illustrate this behavior, we show in Fig. 2 the polynomial reproducing coefficients c i k of a univariate 8th-degree B-spline kernel for i = 0, . . . , 8. It follows that the weight of samples that are away from the origin are considerably larger than the weight of the samples near the origin so that samples near the image borders are dominated noise. In addition, the amplified noise corrupts the numerical moments in a way that increases with the order of the moments. Figure 3 illustrates the numerical reconstruction performance of a typical fourthdegree algebraic shape using our BMR algorithm. We compare our approach against the CMR and RMR methods. In this experiment, white Gaussian noise was added to the image samples with different noise levels using the MATLAB function awgn. As shown in the figure, CMR is very unstable to noise. The comparison shows that our method is significantly more stable than CMR and it performs comparably with RMR.
We remark that, in Fig. 3 and similarly in Figs. 5, 6 and 7, our presentation of the numerical results follows the convention typically adopted by other authors in displaying the reconstruction error. If the panel appears (almost) entirely black, this indicates that the reconstruction error is practically zero; this contrasts with panels where there is a visible error.
To provide a more extensive quantitative comparison of our BMR algorithm with respect to CMR and RMR, we randomly generated 100 bounded algebraic shapes of degree four. Figure 4 compares the reconstruction performance from 47 × 47 noisy samples with SNR = 50 dB. As a performance metric, we used the Dice Similarity Coefficient (DSC) that is commonly used to assess the performance in binary segmentation (images considered here are also binary). This is defined as where T P = true positive pixel, F P = false positive pixel and F N = false negative pixel. DSC ranges between 0 and 1, with DSC = 1 describing perfect reconstruction. The figure shows that CMR is significantly more unstable than our BMR method, while BMR and RMR performs comparably with BMR giving more consistent results overall.
In Fig. 5, we illustrate the sensitivity of our reconstruction approach on the downsampling rate of the original image. As expected, the reconstruction performance Fig. 3 Image reconstruction from noisy samples. Top: algebraic shape of degree 4 (512×512 pixels). Second row: 39 × 39 noisy uniform samples with SNR = 100, 50, and 30 dB. Third-fifth rows: corresponding absolute values of reconstruction error using conventional moment reconstruction (CNR), reconstruction using Bernstein moments (BMR) and regularized reconstruction from [9] (RMR) tends to decrease for higher downsampling rate and this behavior is especially pronounced for the reconstruction algorithm based on the power basis. By contrast, our BMR approach remains stable even for very high rates of downsampling, similarly to RMR.
We remark that, since the RMR algorithm applies an iterative regularized reconstruction strategy, its computational cost is significantly higher than our BMR algorithm. More precisely, it takes our algorithm 22.63 s to process 100 images of size 512×512 as compared to RMR that takes 359.60 s (about 16 times more). Unlike our approach, the computing time of the RMR algorithm may vary from image to image. We also observe that the available numerical code of the RMR is specifically developed to deal with algebraic polynomial of order 4 and is optimized for a SNR Fig. 4 Image reconstruction from noisy samples. Boxplots compare the reconstruction performance of different reconstruction methods on 100 algebraic shapes of degree 4 from 47 × 47 noisy samples (SNR = 50 dB) using the Dice Similarity Coefficient as performance metric = 50 dB. By contrast, our algorithm can be applied to algebraic curves of any degree and with any noise level without any changes.
We also tested our algorithm on the reconstruction of simple shapes that were hand-drawn hence they are not algebraic domains (even though some of them are close to algebraic ones). Results in Fig. 6 show that our BMR algorithm performs well also on these examples. We recall that our algorithm requires to set the degree of algebraic shape to be recovered. Since this parameter is unknown in this case, we set up a simple method to search for the degree of the algebraic shape that best fits the samples as follows. To generate the reconstructions, we ran our algorithm by sequentially choosing degrees = 2, 4, 6, 8 and 10; we next selected the reconstruction with Fig. 5 Sensitivity of reconstruction to sampling rate. Absolute error of the reconstruction from noisy samples with SNR = 50 db at different sampling rates (first column) using reconstruction from conventional moments (CMR), Bernstein moments (BMR) and regularized reconstruction from [9] (RMR) Fig. 6 Reconstruction of non-algebraic shapes. Top: non-algebraic shapes. Bottom: Absolute error of the reconstruction from a 47×47 downsampled version of the image. We assumed polynomial approximations of degree 4, 6 and 6 respectively the highest DSC score. Results on the simple shapes considered in Fig. 6 show that we obtain very satisfactory results using algebraic shapes of relatively low degree. We could not compare our result with RMR in this case since the available code is limited to algebraic curves of degree 4.
We finally examined the sensitivity of our algorithm to the choice of the sampling kernel. Results in Fig. 7 show that, unlike the B-spline kernel, other choices including Fig. 7 Sensitivity of reconstruction to sampling kernel. Algebraic shapes of size of degree four (column 1) were reconstructed using our BMR algorithm from 32 × 32 noisy samples with SNR = 100 dB using different sampling kernels. Absolute error of reconstruction using B-spline, Daubechies and pseudospline kernels is shown in columns 2, 3 an 4 respectively the Daubechies and pseudo-spline kernels are significantly less stable. We attribute this phenomenon to the larger support size of the latter kernels.

Conclusion
We have presented a new strategy for the reconstruction of binary images associated with algebraic shapes using a small number of samples. The main novelty of our approach lies in the use of Bernstein polynomials to represent the image moments. By combining this new representation with the sampling of images using a polynomialgenerating refinable kernel we derive a new reconstruction algorithm that is robust to noise, simple to implement and very fast to compute. Future research includes the investigation of alternative image representations such as the Chebyshev basis which is also known to offer good stability properties.