Sparse image reconstruction on the sphere: implications of a new sampling theorem
Abstract
We study the impact of sampling theorems on the fidelity of sparse image reconstruction on the sphere. We discuss how a reduction in the number of samples required to represent all information content of a bandlimited signal acts to improve the fidelity of sparse image reconstruction, through both the dimensionality and sparsity of signals. To demonstrate this result we consider a simple inpainting problem on the sphere and consider images sparse in the magnitude of their gradient. We develop a framework for total variation (TV) inpainting on the sphere, including fast methods to render the inpainting problem computationally feasible at highresolution. Recently a new sampling theorem on the sphere was developed, reducing the required number of samples by a factor of two for equiangular sampling schemes. Through numerical simulations we verify the enhanced fidelity of sparse image reconstruction due to the more efficient sampling of the sphere provided by the new sampling theorem.
pheres, harmonic analysis, sampling methods, compressive sensing.
1 Introduction
\PARstartImages are observed on a spherical manifold in many fields, from astrophysics (e.g. [1]) and biomedical imaging (e.g. [2]), to computer graphics (e.g. [3]) and beyond. In many of these settings inverse problems arise, where one seeks to recover an unknown image from linear measurements, which may be noisy, incomplete or acquired through a convolution process, for example. Such inverse problems are typically solved by assuming some regularising prior on the unknown image to be recovered.
Sparsity priors have received a lot of attention recently, since (a) they have been shown to be an effective and versatile approach for representing many realworld signals and (b) a sound theoretical foundation is provided by the emerging and rapidly evolving theory of compressive sensing [4, 5, 6]. On the sphere, compressive sensing has been considered for signals sparse in the spherical harmonic domain [7], however a general theoretical framework does not yet exist for signals sparse in spatially localised representations. Nevertheless, sparse image reconstruction on the sphere in alternative representations, such as a set of overcomplete dictionaries, may still be considered; indeed, such an approach has been shown to be very effective [8].
Although compressive sensing goes beyond Nyquist sampling, the Nyquist limit nevertheless defines the benchmark from which compressive sensing improvements are relative. Compressive sensing results are thus tightly coupled to the underlying sampling theorem on the manifold of interest. On the sphere, unlike Euclidean space, the number of samples required in the harmonic and spatial domains differ, with different sampling theorems on the sphere requiring a different number of samples in the spatial domain. Consequently, the sampling theorem adopted influences the performance of sparse signal reconstruction on the sphere. Studying the impact of sampling theorems on the sphere on the performance of sparse signal reconstruction is the focus of the current article.
When considering signal priors that incorporate spatially localised information (for example directly in real space, in the magnitude of the gradient of signals, or through a wavelet basis or overcomplete dictionary), the sampling theorem that is adopted becomes increasingly important. Recently, a new sampling theorem on the sphere was developed by two of the authors of the current article for equiangular sampling schemes [9], reducing Nyquist sampling on the sphere by a factor of two compared to the canonical approach [10, 11]. The reduction in the number of samples required to represent a bandlimited signal on the sphere has important implications for sparse image reconstruction.
To gain some intuition regarding these implications, we appeal to
standard compressive sensing results in Euclidean space, where the
ratio of the number of measurements required to reconstruct a
sparse image, to its dimensionality , goes as [5, 12], where is the sparsity
measure of the image (i.e. the number of nonzero coefficients in some
sparse representation).
Both of these quantities, dimensionality and sparsity, are related to the number of samples required to capture all information content of the underlying signal, as prescribed by the adopted sampling theorem. Spatial dimensionality is given identically by the number of samples of the sampling theorem. For any sparse representation of an image that captures spatially localised information, the sparsity of the signal is also directly related to spatial sampling. For example, in a wavelet representation, wavelets are located on each sample point. A less dense dictionary of wavelet atoms required to span the space will lead to a more sparse representation of images when the sparsity is computed in an analysis approach, i.e. as the number of nonzero projections of the signal onto the wavelet atoms. We concentrate on such analysis priors here, as suggested by the recent evolution of compressive sensing with redundant dictionaries [12]. This argument can be extended to sparsity in the gradient and, in fact, all sparsity measures that capture spatially localised signal content. Consequently, for images sparse in a spatially localised representation, the ability to represent a bandlimited signal on the sphere with fewer samples while still capturing all of its information content will improve the fidelity of sparse image reconstruction by enhancing both the dimensionality and sparsity of signals.
In this article we study the implications of a new sampling theorem [9] for sparse image reconstruction on the sphere. We verify the hypothesis that a more efficient sampling of the sphere, as afforded by the new sampling theorem [9], enhances the fidelity of sparse image reconstruction through both the dimensionality and sparsity of signals. To demonstrate this result we consider a simple inpainting problem, where we recover an image on the sphere from incomplete spatial measurements. We consider images sparse in the magnitude of their gradient, as an illustration of the general setting, and develop a framework for total variation (TV) inpainting on the sphere. Solving these problems is computationally challenging; hence we develop fast methods for this purpose. Our framework is general and is trivially extended to other sparsity priors that incorporate spatially localised information.
The remainder of the article is structured as follows. In Section 2 we concisely review the harmonic structure of the sphere and corresponding sampling theorems. We develop a framework for TV inpainting on the sphere in Section 3. In Section 4 we describe algorithms for solving the optimisation problems on the sphere that arise in our TV inpainting framework. Numerical simulations are performed in Section 5, showing the enhanced fidelity of sparse image reconstruction provided by a more efficient sampling of the sphere. Concluding remarks are made in Section 6.
2 Sampling on the Sphere
A sampling theorem on the sphere states that all information in a (continuous) bandlimited signal is captured in a finite number of samples in the spatial domain. Since a (continuous) bandlimited signal on the sphere may be represented by a finite harmonic expansion, a sampling theorem on the sphere is equivalent to an exact prescription for computing a spherical harmonic transform from a finite set of spatial samples. In this section we review the harmonic structure of the sphere, before discussing sampling theorems on the sphere.
2.1 Harmonic structure of the sphere
We consider the space of square integrable functions on the sphere , with the inner product of defined by
where is the usual invariant measure on the sphere and denote spherical coordinates with colatitude and longitude . Complex conjugation is denoted by the superscript . The canonical basis for the space of square integrable functions on the sphere is given by the spherical harmonics , with natural , integer and . Due to the orthogonality and completeness of the spherical harmonics, any square integrable function on the sphere may be represented by its spherical harmonic expansion
(1) 
where the spherical harmonic coefficients are given by the usual projection onto each basis function:
Throughout, we consider signals on the sphere bandlimited at , that is signals such that , , in which case the summation over in (1) may be truncated to the first terms. Finally, note that the harmonic coefficients of a real function on the sphere satisfy the conjugate symmetry relation , which follows directly from the conjugate symmetry of the spherical harmonics.
2.2 Sampling theorems on the sphere
Sampling theorems on the sphere describe how to sample a bandlimited signal so that all information is contained in a finite number of samples . Moreover, a sampling theorem on the sphere effectively encodes an exact quadrature rule for the integration of bandlimited functions [10, 9]. We denote the concatenated vector of spatial measurements by and the concatenated vector of harmonic coefficients by . The number of spatial and harmonic elements, and respectively, may differ (and in fact do differ for all known sampling theorems on the sphere).
Before discussing different sampling theorems on the sphere, we define a generic notation to describe the harmonic transform corresponding to a given sampling theorem. A sampling theorem describes how to compute the spherical harmonic transform of a signal exactly. Since the spherical harmonic transform and inverse are linear, we represent the forward and inverse transform by the matrix operators and respectively. The spherical harmonic coefficients of a sampled signal (i.e. image) on the sphere are given by the forward transform
while the original signal is recovered from its harmonic coefficients by the inverse transform
Different sampling theorems then differ in the definition of , and the number of spatial samples . By definition, all sampling theorems give exact spherical harmonic transforms, implying , where is the identity matrix. However, for all sampling theorems on the sphere the number of samples required in the spatial domain exceeds the number of coefficients in the harmonic domain (i.e. ), hence . Consequently, for the sample positions of a sampling theorem, an arbitrary set of sample values does not necessarily define a bandlimited signal (contrast this to the discrete Euclidean setting where a finite set of samples uniquely defines a bandlimited signal). Note also that the adjoint inverse (forward) spherical harmonic transform differs to the forward (inverse) spherical harmonic transform in the discrete setting.
For an equiangular sampling of the sphere, the Driscoll & Healy (DH)
[10] sampling theorem has become the standard,
requiring
samples on the sphere to represent exactly a signal bandlimited in
its spherical harmonic decomposition at . Recently, a new
sampling theorem for equiangular sampling schemes has been developed
by McEwen & Wiaux (MW) [9], requiring only
samples to represent a bandlimited signal exactly.
No sampling theorem on the sphere reaches the optimal number of
samples suggested by the dimension of a bandlimited signal
in harmonic space (although the MW sampling theorem comes closest to
this bound). The MW sampling theorem therefore achieves a more
efficient sampling of the sphere, with a reduction by a factor of
approximately two in the number of samples required to represent a
bandlimited signal on the sphere.
Fast algorithms have been developed to compute forward and inverse
spherical harmonic transforms rapidly for both the DH
[10, 11] and MW [9] sampling
theorems. These fast algorithms are implemented, respectively, in the
publicly available SpharmonicKit
3 Sparse Image Reconstruction on the Sphere
A more efficient sampling of a bandlimited signal on the sphere, as afforded by the MW sampling theorem, improves the quality of sparse image reconstruction for images that are sparse in a spatially localised measure. To demonstrate this result we consider a simple inpainting problem on the sphere and consider images sparse in the magnitude of their gradient. We develop a framework for total variation (TV) inpainting on the sphere, which relies on a sampling theorem and its associated quadrature rule to define a discrete TV norm on the sphere. Firstly, we define the discrete TV norm on the sphere, before secondly defining finite difference gradient operators on the sphere. Thirdly, we discuss the TV inpainting problem.
3.1 TV norm on the sphere
We define the discrete TV norm on the sphere by
(2) 
where and index the equiangular samples in and respectively, with the number of samples associated with a given sampling theorem denoted in each dimension by and respectively. The discrete magnitude of the gradient is defined by
(3) 
to approximate the continuous magnitude of the gradient
by finite differences. The finite difference
operators and are defined explicitly
in the following subsection.
The contribution to the TV norm from the magnitude of the gradient for
each pixel value is weighted by the quadrature weights
of the sampling theorem adopted in order to
approximate continuous integration.
3.2 Gradient operators on the sphere
The finite difference operators and defined on the sphere appear in the definition of the discrete magnitude of the gradient given by (3), and thus are required to compute the discrete TV norm on the sphere. Furthermore, as we shall see, to solve the TV inpainting problems outlined in the following subsection, the adjoints of these operators are also required. We define these operators and adjoints explicitly here.
The operator is defined samplewise by
with adjoint
Note that this definition is identical to the typical definition of the corresponding operator on the plane [13]. The operator is defined samplewise by
with adjoint
Since the sphere is periodic in , we define the corresponding finite difference operator to also be periodic. The finite difference operator and adjoint in therefore differ to the typical definition on the plane [13].
The TV norm on the sphere may then be seen as the sum of the magnitude of the weighted gradient
where
for
The weighted gradient operator is defined by
where the weighted finite difference operators are defined by
and
Notice how the inclusion of the weights regularises the term that arises from the definition of the gradient on the sphere, eliminating numerical instabilities that this would otherwise cause. If , corresponding to the South pole of the sphere, then and thus we define to avoid dividing by . Note that the MW sampling theorem includes a sample on the South pole, while the DH sampling theorem does not (neither sampling theorem includes a sample on the North pole). The adjoint weighted gradient operator is then applied as
where the adjoint operators and follow trivially from and .
3.3 TV inpainting on the sphere
We consider the measurement equation
where noisy real measurements of the underlying real image on the sphere are made. The matrix implementing the measurement operator represents a uniformly random masking of the image, with one nonzero, unit value on each row specifying the location of the measured datum. The noise is assumed to be independent and identically distributed (iid) Gaussian noise, with zero mean and variance . We assume that the image is sparse in the norm of its gradient and thus attempt to recover from measurements by solving the following TV inpainting problem directly on the sphere:
(4) 
The square of the residual noise follows a scaled distribution with degrees of freedom, i.e. . Consequently, we choose to correspond to the ()th percentile of this distribution, giving a probability that pure noise produces a residual noise equal to or smaller than the observed residual. Note that the data constraint in (4) is given by the usual norm, which is appropriate for Gaussian noise on a discrete set of measurements. Although we consider bandlimited signals, we have not imposed this constraint when solving (4). Consequently, will not necessarily be bandlimited at and we impose this constraint on the solution by performing a forward and inverse spherical harmonic transform: , where the bandlimiting operator is defined by .
As discussed already, for images sparse in a measure that captures spatially localised information, such as the TV norm, a more efficient sampling of the signal enhances sparsity. Furthermore, when recovering signals in the spatial domain directly, the dimensionality of the signal is also enhanced by a more efficient sampling. These two effects both act to improve the fidelity of sparse image reconstruction. Thus, the more efficient sampling of the MW sampling theorem when compared to the DH sampling theorem will improve the fidelity of sparse image reconstruction when solving the TV inpainting problem given by (4). We verify these claims with numerical experiments in Section 5.
No sampling theorem on the sphere reaches the optimal number of samples in the spatial domain suggested by the dimensionality of the signal in the harmonic domain. We may therefore optimise the dimensionality of the signal that we attempt to recover by recovering its harmonic coefficients directly. We do so by solving the following TV inpainting problem in harmonic space:
(5) 
We impose reality of the recovered signal by explicitly imposing conjugate symmetry in harmonic space through the conjugate symmetry extension operator , where . The full set of harmonic coefficients of are given by , where are the harmonic coefficients for the spherical harmonic azimuthal index nonnegative only. The image on the sphere is then recovered from its harmonic coefficients by . By solving the TV inpainting problem directly in harmonic space, we naturally recover a signal bandlimited at .
When solving the TV inpainting problem (5) directly in harmonic space, the dimensionality of the recovered signal is optimal and identical for both sampling theorems. However, the sparsity of the signal with respect to the TV norm remains enhanced for the MW sampling theorem when compared to the DH sampling theorem. Consequently, the MW sampling theorem will improve the fidelity of sparse image reconstruction when solving the TV inpainting problem given by (5), although through sparsity only and not also dimensionality. We verify these claims with numerical experiments in Section 5.
Note that if a bandlimit constraint were explicitly imposed in problem (4), then the two problems would be equivalent, however, this would involve applying the bandlimiting operator , complicating the problem and increasing the computational cost of finding a solution, while providing no improvement over (5). In the current formulation of these two optimisation problems, problem (4) has the advantage of simplicity, while problem (5) is the simplest formulation that optimises dimensionality.
4 Algorithms
We solve the TV inpainting problems on the sphere given by (4) and (5) using iterative convex optimisation methods. Solving the TV inpainting problem in harmonic space poses two challenges as we go to high bandlimits (i.e. highresolution). Firstly, we require as an input to the convex optimisation algorithm an upper bound on the inverse transform operator norm, which is challenging to compute at highresolution. We describe a method to compute the operator norm at highresolution, which, crucially, does not require an explicit computation of . Secondly, the inverse spherical harmonic transform and its adjoint operator must be applied repeatedly in the iterative algorithm. Fast algorithms are essential to perform forward and inverse spherical harmonic transforms at highresolution and have been developed for both the DH [10, 11] and MW [9] sampling theorems. To solve the inpainting problem at highresolution we also require a fast adjoint inverse transform. We thus develop fast algorithms to perform the adjoint forward and adjoint inverse spherical harmonic transforms corresponding to the MW sampling theorem. Since we predict the MW sampling theorem to be superior to the DH sampling theorem for sparse image reconstruction on the sphere (a prediction that is validated by numerical experiments performed at low bandlimits (i.e. lowresolution) in Section 5), we develop fast adjoint algorithms for the MW sampling theorem only. These methods then render the computation of solutions to the TV inpainting problems feasible at highresolution for the MW sampling theorem.
4.1 Convex optimisation
We apply the DouglasRachford proximal splitting algorithm
[14] to solve the convex optimisation problems
(4) and
(5).
The DouglasRachford algorithm [14] is based on a splitting approach that requires the computation of two proximity operators [15]. In our case, one proximity operator is based on the TV norm and the other on the data constraint .
In the case of an image on the plane, the proximity operator based on the TV norm may be computed using, for example, the method described in [13] or in [16]. For an image on the sphere, the same methods can be used after introducing the following modifications. In [16] the algorithm to compute the proximity operator of the TV norm is described in terms of a linear operator , its adjoint , and two projections onto a set and a set . In our case, the linear operator and its adjoint may be redefined as
and
where the set is the set of weighted gradientpairs such that and is simply given by the space of the recovered vector .
The second proximity operator, related to the data constraint , is computed using the method described in [17] directly.
4.2 Operator norm bound
The convex optimisation algorithm requires as input upper bounds for the norms of the operators that appear in the problem. The calculation of these norms is in most cases straightforward, however the calculation of the inverse spherical harmonic transform operator norm, defined by
can prove problematic. At lowresolution may be computed explicitly, however this is not feasible at highresolution since even computing and storing explicitly is challenging.
We develop a method here to estimate this norm for the MW sampling theorem without computing explicitly. We seek a sampled function on the sphere that maximises , while satisfying the constraint . By the Parseval relation and the sampling theorem on the sphere, this constraint may be rewritten:
where contains samples of , sampled at a resolution sufficient to represent , i.e. corresponding to bandlimit (so that an exact quadrature may be used to evaluate from a discrete set of samples), is the matrix with corresponding quadrature weights along its diagonal, and where . Since we know that the quadrature weights for the MW sampling theorem are closely approximated by [9], the signal that maximises while satisfying the constraint has its energy centred as much as possible on the South pole since this is where the quadrature weights are smallest (recall that the MW sampling scheme does not contain a sample on the North pole). This signal is given by the bandlimited Dirac delta function centred on the South pole (see e.g. [18] for the definition of the bandlimited Dirac delta function on the sphere). The spherical harmonic coefficients of this bandlimited Dirac delta function are given by
where is a normalisation factor chosen to ensure and is the Kronecker delta symbol. The norm of the inverse spherical harmonic transform operator may then be computed by , which, crucially, does not require an explicit computation of , merely its application.
In Figure 1 we compute by the method outlined here and from explicitly, for lowresolution. We find that the method to estimate the norm of the inverse spherical harmonic transform operator outlined here estimates the actual norm very well.
We also derived an upper bound for the norm of this operator for the MW sampling theorem. However, the bound we derived is not tight and we found empirically that the method outlined here to estimate the norm itself, rather than a bound, is very accurate and improved the performance of the optimisation algorithm considerably when compared to a nontight bound. Although we do not prove so explicitly, we conjecture that the method outlined here gives the inverse transform operator norm exactly.
4.3 Fast adjoint spherical harmonic transforms
Standard convex optimisation methods require not only the application of the operators that appear in the optimisation problem but often also their adjoints. Moreover, these methods are typically iterative, necessitating repeated application of each operator and its adjoint. Thus, to solve optimisation problems that incorporate harmonic transform operators, like the harmonic space TV inpainting problem given by (5), fast algorithms to apply both the operator and its adjoint are required to render highresolution problems computationally feasible.
Here we develop fast algorithms to perform adjoint forward and adjoint
inverse spherical harmonic transforms for the MW sampling theorem.
Although we only require the adjoint inverse transform in this
article, for the sake of completeness we also derive a fast adjoint
forward transform. Similarly, although we only consider scalar
functions in this article, for the sake of completeness we derive fast
adjoint algorithms for the spin setting. A spin function on the sphere
transforms as
under a local rotation by , where the prime denotes
the rotated function. It is important to note that the rotation
considered here is not a global rotation on the sphere but
rather a rotation by in the tangent plane at (see e.g. [9] for further details). In the expressions for the
fast algorithms derived below, the standard scalar case follows simply
by setting . These fast adjoint algorithms are implemented
in the publicly available SSHT
The fast adjoint inverse spherical harmonic transform for the MW sampling theorem follows by taking the adjoint of each stage of the fast inverse transform [9] and applying these in reverse order. The final stage of the fast inverse transform involves discarding outofdomain samples and has adjoint
The second stage of the fast adjoint inverse transform is given by
which may be computed rapidly using fast Fourier transforms (FFTs). The final stage of the fast adjoint inverse transform is given by
where are the Wigner functions evaluated for argument (see e.g. [19]). This final calculation dominates the overall asymptotic complexity of the fast adjoint inverse transform, resulting in an algorithm with complexity .
The fast adjoint forward spherical harmonic transform for the MW sampling theorem follows by taking the adjoint of each stage of the fast forward transform [9] and applying these in reverse order. The first stage of the fast adjoint forward transform is given by
The next stage is given by the (reflected) convolution
which is selfadjoint, followed by the inverse Fourier transform in
which may be computed rapidly using FFTs. The next stage consists of the adjoint of the periodic extension of a function on the sphere performed in the forward transform and is given by
The final stage consists of the Fourier transform in
which may be computed rapidly using FFTs. The first calculation dominates the overall asymptotic complexity of the fast adjoint forward transform, resulting in an algorithm with complexity .
5 Simulations
We perform numerical experiments to examine the impact of a more efficient sampling of the sphere when solving the TV inpainting problems defined in Section 3. Firstly, we perform a lowresolution comparison of reconstruction fidelity when adopting the DH and MW sampling theorems, where the predicted improvements in reconstruction fidelity provided by the MW sampling theorem are verified in practice. Secondly, we perform a single simulation to illustrate TV inpainting at highresolution on a realistic test image.
5.1 Lowresolution comparison on bandlimited images
A test image is constructed from Earth topography data. The original
Earth topography data are taken from the Earth Gravitational Model
(EGM2008) publicly released by the U.S. National
GeospatialIntelligence Agency (NGA) EGM Development
Team.
Measurements of the test image are taken at uniformly random locations on the sphere, as described by the measurement operator , in the presence of Gaussian iid noise with standard deviation . Reconstructed images on the sphere are recovered by solving the inpainting problems in the spatial and harmonic domains, through (4) and (5) respectively, using both the DH and MW sampling theorems, giving four reconstruction techniques. The bound is determined from . We consider the measurement ratios (recall that is the dimensionality of the signal in harmonic space). The measurement ratio corresponds to complete coverage for the MW sampling theorem, i.e. Nyquist rate sampling on the MW grid.
Typical reconstructed images are shown in Figure 3 for the four reconstruction techniques. For each reconstruction technique and measurement ratio , we perform ten simulations for random measurement operators and noise. To quantify the error of reconstruction, we compute the signaltonoiseratio (defined in harmonic space to avoid differences due to the number of samples of each sampling theorem). Note that the standard norm is used in the definition of the SNR given the discrete nature of harmonic space on the sphere. Reconstruction performance, averaged over these ten simulations, is shown in Figure 4.
When solving the inpainting problem in the spatial domain through (4) we see a large improvement in reconstruction quality for the MW sampling theorem when compared to the DH sampling theorem. This is due to the enhancement in both dimensionality and sparsity afforded by the MW sampling theorem in this setting. When solving the inpainting problem in the harmonic domain through (5) we see a considerable improvement in reconstruction quality for each sampling theorem, since we optimise the dimensionality of the recovered signal by going to harmonic space. For harmonic reconstructions, the MW sampling theorem remains superior to the DH sampling theorem due to the enhancement in sparsity (but not dimensionality) that it affords in this setting. All of the predictions made in Section 3 are thus exhibited in the numerical experiments performed in this section. In all cases, the superior performance of the MW sampling theorem is clear.
5.2 Highresolution illustration on a realistic image
In this section we perform a single simulation to illustrate TV inpainting at high resolution. Furthermore, we also consider a more realistic test image. Since we develop fast adjoint algorithms for the MW sampling theorem only (due to its superiority), we therefore use only the MW sampling theorem for the highresolution inpainting simulation performed here.
A highresolution test image is constructed from the same Earth topography data described in Section 5.1. Since the original data are defined in harmonic space, we first simply truncate the harmonic coefficients to yield an image bandlimited at . In practice, acquired images may not necessarily be bandlimited. We thus process the data to construct a test image that is not bandlimited. We construct such a test image defined by height above sealevel (i.e. we threshold all oceans and trenches, while the continents and mountains remain unaltered). The abrupt transition between the oceans and the continents results in an image that is indeed not bandlimited (as verified numerically). Furthermore, the continents and mountainous regions result in a test image that is not highly sparse in its gradient. The resulting realistic test image is shown in Figure 5 (a).
The same measurement procedure as outlined previously is applied to take noisy, incomplete measurements of the data for a range of measurement ratios . The inpainted images are recovered by solving the inpainting problem in harmonic space through (5) using the MW sampling theorem. To solve the inpainting problem for these highresolution simulations we use the estimator of the inverse transform norm described in Section 4.2 and the fast adjoint harmonic transform algorithms defined in Section 4.3. Using these fast algorithms, combined with recent optimisations of the SSHT package, it takes approximately 10 minutes to solve the inpainting problem in harmonic space at on a standard laptop (with a 1.8 GHz Intel Core i7 processor and 4 GB of RAM).
The inpainted images are shown in Figure 5. Since the original realistic test image is not bandlimited, the previous SNR measure (which is defined in harmonic space to avoid a dependence on the number of samples of each sampling theorem) is not a meaningful error metric (computation of the harmonic transform would indeed be affected by uncontrolled aliasing). Instead, we use the analogous SNR measure defined in image space on the sphere, given by , where we recall is the matrix with quadrature weights on its diagonal. Just as for the definition of the TV norm, the inclusion of the weights for signals that are not bandlimited provides only an intuitive approximation to continuous integration. The values for each inpainted image are displayed in Figure 5, from which it is apparent that increases with increasing measurement ratio. Moreover, it is clearly apparent by eye that our TV inpainting framework is effective when applied to realistic images that are not highly sparse in their gradient and that are not bandlimited.






6 Conclusions
The MW sampling theorem, developed only recently, achieves a more efficient sampling of the sphere than the standard DH sampling theorem: without any loss to the information content of the sampled signal, the MW sampling theorem reduces the number of samples required to represent a bandlimited signal by a factor of two for an equiangular sampling. For signals sparse in a spatially localised measure, such as in a wavelet basis, overcomplete dictionary, or in the magnitude of their gradient, for example, a more efficient sampling enhances the fidelity of sparse image reconstruction through both dimensionality and sparsity. When a signal is recovered directly in the spatial domain, the MW sampling theorem provides enhancements in both dimensionality and sparsity when compared to the DH sampling theorem. By recovering the signal directly in harmonic space it is possible to optimise its dimensionality, in which case the MW sampling theorem still provides an enhancement in sparsity but not in dimensionality.
We verified these statements through a simple inpainting problem on the sphere, where we considered images sparse in their gradient. We built a framework and fast methods for total variation (TV) inpainting on the sphere. Using this framework we performed numerical experiments which confirmed our predictions: in all cases, the more efficient sampling provided by the MW sampling theorem improved the fidelity of sparse image reconstruction on the sphere.
Jason McEwen received a B.E. (Hons) degree in Electrical and Computer Engineering from the University of Canterbury, New Zealand, in 2002 and a Ph.D. degree in Astrophysics from the University of Cambridge in 2006.
He held a Research Fellowship at Clare College, Cambridge, from 2007 to 2008, worked as a Quantitative Analyst from 2008 to 2010, and held a position as a Postdoctoral Researcher at Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland, from 2010 to 2011. From 2011 to 2012 he held a Leverhulme Trust Early Career Fellowship at University College London (UCL), where he remains as a Newton International Fellow, supported by the Royal Society and the British Academy. His research interests are focused on spherical signal processing, including sampling theorems and wavelets on the sphere, compressed sensing and Bayesian statistics, and applications of these theories to cosmology and radio interferometry.
Gilles Puy was born in Nevers, France, on August 8, 1985. He received the Engineering degree from the Ecole Supérieure d’Electricité (Supélec), GifsurYvette, France, and the M. Sc. degree in Electrical and Electronics Engineering from the Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland, in 2009.
In 2009, he joined the Signal Processing Laboratory and the Laboratory of functional and metabolic imaging at the EPFL as a Doctoral Assistant, where he is currently working towards the Ph.D. degree in Electrical Engineering. His main research interests include inverse problems, compressed sensing, biomedical imaging and radio interferometry.
JeanPhilippe Thiran received the Electrical Engineering and Ph.D. degrees from the Université catholique de Louvain (UCL), LouvainlaNeuve, Belgium, in 1993 and 1997, respectively.
Since January 2004, he has been an Assistant Professor, responsible for the Image Analysis Group at the Swiss Federal Institute of Technology (EPFL), Lausanne, Switzerland. His current scientific interests include image segmentation, prior knowledge integration in image analysis, partial differential equations and variational methods in image analysis, multimodal signal processing, medical image analysis, including multimodal image registration, segmentation, computerassisted surgery, and diffusion MRI.
Dr. Thiran was CoEditorinChief of Signal Processing (published by Elsevier Science) from 2001 to 2005. He is currently an Associate Editor of the International Journal of Image and Video Processing (published by Hindawi), and member of the Editorial Board of Signal, Image and Video Processing (published by Springer). He was the General Chairman of the 2008 European Signal Processing Conference (EUSIPCO 2008). He is a senior member of the IEEE, and a member of the MLSP and IVMSP technical committees of the IEEE Signal Processing Society.
Pierre Vandergheynst received the M.S. degree in Physics and the Ph.D. degree in Mathematical Physics from the Université catholique de Louvain, LouvainlaNeuve, Belgium, in 1995 and 1998, respectively.
From 1998 to 2001, he was a Postdoctoral Researcher and an Assistant Professor with the Signal Processing Laboratory, Swiss Federal Institute of Technology (EPFL), Lausanne, Switzerland. He is now an Associate Professor at EPFL, where his research focuses on harmonic analysis, sparse approximations, and mathematical image processing with applications to higher dimensional, complex data processing.
Dr. Vandergheynst was coEditorinChief of Signal Processing from 2002 to 2006 and has been an Associate Editor of the IEEE TRANSACTIONS ON SIGNAL PROCESSING since 2007. He has been on the Technical Committee of various conferences and was CoGeneral Chairman of the EUSIPCO 2008 conference. He is the author or coauthor of more than 50 journal papers, one monograph and several book chapters. He is a laureate of the Apple ARTS award and holds seven patents.
Dimitri Van De Ville received the M.S. degree in Engineering and Computer Sciences from Ghent University, Belgium, in 1998, as well as the Ph.D. degree in Computer Science Engineering, in 2002.
He obtained a grant as Research Assistant with the Fund for Scientific Research Flanders Belgium (FWO). In 2002, he joined Prof. M. Unser’s Biomedical Imaging Group at the Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland. In December 2005, he became responsible for the Signal Processing Unit at the University Hospital of Geneva, Geneva, Switzerland, as part of the Centre d’Imagerie Biomédicale (CIBM). He was recently awarded a SNSF professorship from the Swiss National Science Foundation and he currently holds a joint position at the University of Geneva and the EPFL. His research interests include wavelets, sparsity, pattern recognition, and their applications in biomedical imaging, such as functional magnetic resonance imaging.
Dr. Van De Ville served as an Associate Editor for the IEEE Transactions on Image Processing (20062009) and the IEEE Signal Processing Letters (20042006). He is a member of the Bio Imaging and Signal Processing (BISP) TC of the IEEE SPS. Since 2003, he has also been an Editor and Webmaster of The Wavelet Digest. He is cochair of the Wavelets series conferences (2007, 2009), together with V. Goyal and M. Papadakis.
Yves Wiaux received the M.S. degree in Physics and the Ph.D. degree in Theoretical Physics from the Université catholique de Louvain (UCL), LouvainlaNeuve, Belgium, in 1999 and 2002, respectively.
He was a Postdoctoral Researcher at the Signal Processing Laboratories of the Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland, from 2003 to 2008. He was also a Postdoctoral Researcher of the Belgian National Science Foundation (F.R.S.FNRS) at the Physics Department of UCL from 2005 to 2009. He is now a Maître Assistant of the University of Geneva (UniGE), Switzerland, with joint affiliation between the Institute of Electrical Engineering and the Institute of Bioengineering of EPFL, and the Department of Radiology and Medical Informatics of UniGE. His research lies at the intersection between complex data processing (including development on wavelets and compressed sensing) and applications in astrophysics (notably in cosmology and radio astronomy) and in biomedical sciences (notably in MRI and diffusion MRI).
Footnotes
 Typically the mutual coherence of the measurement and sparsifying operators also plays a role [5]. However, in Euclidean space, as on the sphere, discrete inner products can be related to the (unique) continuous inner product (via a sampling theorem). Consequently, the measure of coherence is invariant to the choice of sampling theorem (the coherence is defined through the continuous inner product, which is the same for all sampling theorems). For the purpose of comparing sampling theorems on the sphere, we can thus safely neglect the impact of coherence.
 GaussLegendre (GL) quadrature can also be used to construct an efficient sampling theorem on the sphere, with samples (see e.g. [9]). The MW sampling theorem nevertheless remains more efficient, especially at low bandlimits. Furthermore, it is not so straightforward to define the TV norm on the GL grid since it is not equiangular. Finally, algorithms implementing the GL sampling theorem have been shown to be limited to lower bandlimits and less accurate than the algorithms implementing the MW sampling theorem [9]. Thus, we focus on equiangular sampling theorems only in this article.
 http://www.cs.dartmouth.edu/~geelong/sphere/
 http://www.jasonmcewen.org/
 If the bandlimiting operator were applied to in (2), then the finite summation of (2) would give an exact quadrature for the integral of the continuous function underlying the associated samples of the bandlimited . However, introducing the operator in (2) would make solving the optimisation problems defined subsequently problematic and would also prohibit passing the quadrature weights inside the gradient to eliminate numerical instabilities due to the term. Consequently, we adopt the definition of the discrete TV norm on the sphere given by (2). In any case, numerical experiments have shown that is identical to to machine precision for the particular test images considered in Section 5.1. Thus, the discrete TV norm defined by (2) can be though of as an accurate proxy for .
 We use DouglasRachford splitting since this does not require differentiability of the objective function and allows us to solve constrained optimisation problems where we adopt an indicator function to represent the measurement constraint.
 http://www.jasonmcewen.org/
 These data were downloaded and extracted using the tools available from Frederik Simons’ webpage: http://www.frederik.net.
References
 N. Jarosik, C. L. Bennett, J. Dunkley, B. Gold, M. R. Greason, M. Halpern, R. S. Hill, G. Hinshaw, A. Kogut, E. Komatsu, D. Larson, M. Limon, S. S. Meyer, M. R. Nolta, N. Odegard, L. Page, K. M. Smith, D. N. Spergel, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, “Sevenyear Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Sky Maps, Systematic Errors, and Basic Results,” Astrophys. J. Supp., vol. 192, p. 14, Feb. 2011.
 H. JohansenBerg and T. E. J. Behrens, Diffusion MRI: From quantitative measurement to invivo neuroanatomy. San Diego: Academic Press, 2009.
 R. Ramamoorthi and P. Hanrahan, “A signal processing framework for reflection,” ACM Transactions on Graphics, vol. 23, no. 4, pp. 1004–1042, 2004.
 E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory, vol. 52, no. 2, pp. 489–509, feb 2006.
 E. Candès, “Compressive sampling,” in Proc. Int. Congress Math., ser. Euro. Math. Soc., vol. III, 2006, p. 1433.
 D. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory, vol. 52, no. 4, pp. 1289–1306, apr 2006.
 H. Rauhut and R. Ward, “Sparse recovery for spherical harmonic expansions,” in Proc. SampTA, 2011.
 P. Abrial, Y. Moudden, J.L. Starck, B. Afeyan, J. Bobin, J. Fadili, and M. K. Nguyen, “Morphological component analysis and inpainting on the sphere: application in physics and astrophysics,” J. Fourier Anal. and Appl., vol. 14, no. 6, pp. 729–748, 2007.
 J. D. McEwen and Y. Wiaux, “A novel sampling theorem on the sphere,” IEEE Trans. Sig. Proc., vol. 59, no. 12, pp. 5876–5887, 2011.
 J. R. Driscoll and D. M. J. Healy, “Computing Fourier transforms and convolutions on the sphere,” Advances in Applied Mathematics, vol. 15, pp. 202–250, 1994.
 D. M. J. Healy, D. Rockmore, P. J. Kostelec, and S. S. B. Moore, “FFTs for the 2sphere – improvements and variations,” J. Fourier Anal. and Appl., vol. 9, no. 4, pp. 341–385, 2003.
 E. J. Candès, Y. Eldar, D. Needell, and P. Randall, “Compressed sensing with coherent and redundant dictionaries,” Applied Comput. Harm. Anal., vol. 31, no. 1, pp. 59–73, 2010.
 A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis., vol. 20, no. 1–2, pp. 89–97, 2004.
 P. Combettes and J.C. Pesquet, “A douglasrachford splitting approach to nonsmooth convex variational signal recovery,” IEEE J. Selected Top. in Sig. Proc., vol. 1, no. 4, pp. 564–574, dec 2007.
 P. Combettes and J.C. Pesquet, Proximal splitting methods in signal Pprocessing. New York: Springer, 2011.
 A. Beck and M. Teboulle, “Fast gradientbased algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. Image Proc., vol. 18, no. 11, pp. 2419–2434, 2009.
 M. J. Fadili and J. L. Starck, “Monotone operator splitting for optimization problems in sparse recovery,” in ICIP, nov 2009, pp. 1461–1464.
 F. J. Simons, F. A. Dahlen, and M. A. Wieczorek, “Spatiospectral concentration on a sphere,” SIAM Rev., vol. 48, no. 3, pp. 504–536, 2006.
 D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, Quantum theory of angular momentum. Singapore: World Scientific, 1989.